Journal  of  Power  Sources  196  (2011)  6534-6553 


ELSEVIER 


Contents  lists  available  at  ScienceDirect 

Journal  of  Power  Sources 

journal  homepage:  www.elsevier.com/locate/jpowsour 


in 

SllbudtS 


Electrochemical  modeling  of  single  particle  intercalation  battery  materials  with 
different  thermodynamics 

Wei  Lai* 

Department  of  Chemical  Engineering  and  Materials  Science,  Michigan  State  University,  East  Lansing,  MI  48824,  USA 


ARTICLE  INFO 


ABSTRACT 


Article  history: 

Received  29  January  201 1 

Received  in  revised  form  21  March  2011 

Accepted  26  March  2011 

Available  online  6  April  2011 


Keywords: 

Single  particle 
Intercalation  battery 
Regular  solution 
Phase  transformation 
Poisson-Nernst-Planck 
Current-voltage  relation 


Current-voltage  relations  of  a  single  intercalation  battery  electrode  particle  were  modeled  with  the 
step  current,  step  voltage,  linear  sweep  voltage,  and  sinusoidal  current  signals.  A  solid  solution  with 
constant  diffusivity,  a  solid  solution  with  variable  diffusivity,  and  a  phase  transformation  material  were 
considered  for  their  thermodynamic  and  kinetic  evaluations  based  on  the  regular  solution  model  and  the 
generalized  Poisson-Nernst-Planck  equations.  The  numerical  simulation  results  were  compared  with 
known,  small-signal  solutions  and  experimental  data  throughout  the  article. 
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1.  Introduction 

Electrode  particles  in  intercalation  batteries  work  by  the 
insertion/intercalation  and  extraction/deintercalation  of  ions  and 
electrons  to  and  from  the  host.  In  the  case  of  lithium  intercalation, 
lithium  ions  enter  the  interstitial  polyhedral  sites  while  electrons 
bind  to  the  transition  metal  sites  as  polarons  or  band  electrons.  Not 
only  lithium  ions  and  electrons  occupy  different  sites  in  the  host, 
they  also  come  from  different  external  sources,  e.g.  electrolyte  for 
lithium  ions  and  current  collector  for  electrons. 

The  ambipolar  charge  transport  of  ions  and  electrons  can  be 
modeled  by  the  Poisson-Nernst-Planck  (PNP)  equations  composed 
of  the  Nernst-Planck,  continuity  equations,  and  Poisson  equation. 
The  PNP  equations  have  been  widely  used  in  liquid  electrochem¬ 
istry,  semiconductors,  and  ionic  channels  in  biological  membranes 
[1-4],  among  others.  Another  popular  view  of  lithium  ion  and  elec¬ 
tron  insertion  is  to  treat  it  as  the  equivalent  of  the  insertion  of  a 
neutral  lithium  species,  which  has  been  modeled  by  Fields  laws  [5]. 
Both  the  classical  PNP  equations  and  Fields  laws  are  partial  differen¬ 
tial  equations  that  track  concentrations  and  mass  fluxes.  This  makes 
them  inconvenient  for  electrochemical  modeling,  which  typically 
concerns  the  current-voltage  relationship.  Classical  PNP  equations 
can  be  reformulated  as  equations,  that  involve  current  (density) 
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and  voltage  variables  [6-13],  which  also  have  the  apparent  form  of 
Fields  laws  [11,12]. 

Electrochemical  modeling  focuses  on  current-voltage  rela¬ 
tions,  where  common  electrical  signals  consist  of  step  currents, 
step  voltages,  linear  sweep/scanning  voltages,  and  sinusoidal 
current/voltages.  All  these  signals  can  be  found  in  classical  elec¬ 
trochemistry  books  by  Bard  et  al.  [14,15].  Step  current  or  voltage 
techniques  are  called  constant-current  chronopotentiometry  and 
potential  step  chronoamperometry  (PSCA),  respectively,  or  simply 
chronopotentiometry  or  chronoamperometry.  If  a  single  current 
or  voltage  step  is  applied,  the  technique  is  commonly  known 
as  constant-current  or  constant-voltage  charging/discharging  in 
the  battery  field  [16].  If  the  steps  are  performed  sequentially, 
chronopotentiometry  and  chronoamperometry  are  classified  as 
galvanostatic  intermittent  titration  technique  (GITT)  [17]  and 
potentiostatic  intermittent  titration  technique  (PITT)  [18],  respec¬ 
tively,  by  Huggins  et  al.  The  technique  with  a  linear  sweep  voltage 
is  called  linear  sweep  voltammetry  (LSV).  If  the  sweeping  is 
performed,  both  anodically  and  cathodically,  it  is  called  cyclic 
voltammetry  (CV).  Finally,  the  technique  involving  sinusoidal  elec¬ 
trical  signal  is  usually  called  impedance  spectroscopy  (IS)  [19], 
dielectric  spectroscopy  (DS)  [20],  or  electrochemical  impedance 
spectroscopy  (EIS)  when  applied  to  electrochemical  systems  [21  ]. 

While  it  is  experimentally  routine  to  apply  several  of  the 
above-mentioned  techniques  simultaneously  for  the  investigation 
of  battery  materials,  (e.g.  in  Ref.  [22]),  some  studies  on  numerical 
simulation  only  focus  on  one  technique,  such  as  constant  current 
[23-25],  potential  step  [26,27],  linear  sweep  voltammetry  [28,29], 
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Nomenclature 

b  mobility 

c  concentration 

c°  reference  or  maximum  concentration 

Qchem  chemical  capacitance 

D  chemical  diffusivity 

e  electron  charge 

/  dimensionless  frequency 

/exp  experimental  frequency 

g  dimensionless  interaction  parameter 

j  dimensionless  current  density 

J  current  density 

f ns  mass  flux 

1<b  Boltzmann  constant 

ks  jump  parameter  in  the  solid 

ki  jump  parameter  in  the  liquid 

m  planar,  cylindrical,  and  spherical  symmetry  param¬ 

eter 

r  length 

R  resistance 

R[nt  interfacial  resistance 

t  time 

t0  characteristic  time 

T  absolute  temperature 

V  voltage 

x  position 

X  dimensionless  concentration 

z  charge  number 

Greek  letters 

a  apparent  symmetry  factor 

8Sl  solid  |  liquid  interfacial  length 

£0  vacuum  permittivity 

sr  relative  permittivity 

0  electrical  potential 

y  dimensionless  interfacial  energy 

/x  chemical  potential 

/x*  reduced  chemical  potential 

/x  electrochemical  potential 

/x*  reduced  electrochemical  potential 

o  conductivity 

r  dimensionless  time 

co  angular  frequency 

0  dimensionless  voltage 


and  impedance  spectroscopy  [30].  In  these  studies  Fields  laws  are 
almost  invariably  used  with  constant  diffusivity.  The  application  of 
constant  diffusivity  is  based  on  the  assumption  of  “small-signal” 
perturbations;  however,  it  is  not  clear  what  magnitude  should 
be  considered  “small”.  In  addition,  kinetic  properties  (diffusivity) 
and  thermodynamic  properties  (chemical  potential  or  voltage)  are 
usually  treated  as  two  separate  sets  of  data,  while  it  is  recog¬ 
nized  that  the  evolution  of  systems  (kinetics)  is  determined  by 
the  thermodynamic  information  under  the  context  of  irreversible 
thermodynamics  [31]. 

This  work  seeks  to  provide  a  comprehensive  numerical  study  on 
the  current-voltage  relationship  of  a  single  intercalation  electrode 
particle  with  different  thermodynamic  properties.  The  thermo¬ 
dynamic  model  utilized  a  regular  solution  model  or  Frumkin 
adsorption  isotherm,  which  will  be  discussed  in  Section  2.2.  This 
study  will  consider  three  different  cases,  including  a  solid  solution 
with  constant  diffusivity,  a  solid  solution  with  variable  diffusivity, 
and  a  phase  transformation  material.  The  evolution  of  the  system,  in 


the  form  of  the  transport  of  ions  and  electrons,  will  be  described  by 
the  generalized  PNP  equations  in  Section  2.1.  Boundary  conditions 
and  dimensionless  forms  of  involved  partial  different  equations  are 
described  in  Sections  2.3  and  2.4.  Numerical  simulation  results  for 
the  four  common  electrical  signals  such  asthe  step  current,  step 
voltage,  linear  sweep  voltage,  and  sinusoidal  current,  are  presented 
in  Section  3. 

2.  Models 

The  detailed  discussion  of  the  generalized  PNP  equations  and 
the  regular  solution  model  can  be  found  in  references  [11,12]  and 
is  included  here  for  completeness. 

2.1.  Generalized  Poisson-Nernst-Planck  (PNP)  equations 

The  classical  PNP  equations  [1-4]  are  usually  formulated  as 
coupled  partial  differential  equations  that  track  concentration  q, 
electrical  potential  0,  and  mass  flux J™s  as  follows: 

r = (i) 

£  +  V -J™  =  0,  (2) 

-ere0A<p  =  y>eq,  (3) 

i 

where  Dz  is  the  chemical  diffusivity,  z,-  is  the  charge  number,  e  is 
the  elementary  electron  charge,  and  er  and  s0  are  relative  and  vac¬ 
uum  permittivity,  respectively.  The  conductivity  <7,  is  the  product 
of  concentration  q  and  mobility  b[  as  follows: 

a,-  =  (Zjefqbj.  (4) 

For  the  dilute  solution,  these  equations  are  known  as  drift-diffusion 
equations  in  semiconductor  devices  [3]. 

These  constitutive  equations  in  electrochemistry  can  be  refor¬ 
mulated  as  generalized  Poisson-Nernst-Planck  equations  [11,12] 


Ji  =  -0i-v  A*, 

(5) 

3zx* 

cihem-gf- +v  Ji  =  °> 

(6) 

9 

Jdis  =  -Yt{£reoV4>), 

(7) 

Jt  =Jdis  T  'y  Ji' 

(8) 

This  formulation  was  primarily  discussed  under  the  context  of 
dilute  solutions  and  impedance  spectra  reported  in  the  past  [6-10]. 

The  correlation  between  the  two  formulations  is  given  by  sev¬ 
eral  physically  intuitive  but  less  commonly  used  parameters.  First, 
the  current  density  Jit  reduced  electrochemical  potential  /x?,  and 
reduced  chemical  potential  /x*  are  used  in  Eqs.  (5)-(8),  instead  of 
the  conventional  mass  flux  J™s,  electrochemical  potential  /Xj,  and 
chemical  potential  /xz  found  in  Eqs.  (l)-(3).  The  relation  is  given  by 
the  following: 

Ji  =  zieJ'/1S<  Ai  =  zieA:  >  Mi  =  Zfill*.  (9) 

The  electrochemical  potential  /q  is  the  sum  of  chemical  potential 
/xz  and  electrical  energy  Z;e0  as  follows: 

A,-  =  +  Z,e0.  (10) 

Thus,  the  reduced  electrochemical  potential  has  the  form: 

Ai=/t*  +  </>-  (ii) 
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Fig.  1.  Equivalent  circuit  for  a  two-charge  carrier  (positive  +  and  negative  -)  particle  under  electroneutrality  conditions.  0  and  r0  are  at  the  center  and  surface  of  the  particle, 
respectively.  Differential  resistor  ( dR+  and  dR_)  and  capacitor  ( dC £hem)  elements  are  shown  along  with  surface  voltages  A+(ro)  and  /Z*  (r0). 


Second,  volumetric  chemical  capacitance  c^hem  represents  the 
change  of  volumetric  electrical  charge  z2ec2-  upon  a  change  of 
reduced  chemical  potential  /x*: 


Qchem 


z.e 


d/x*' 


(12) 


Third,  chemical  diffusivity  D2  is  defined  as  the  following: 


Dj  = 


■* chem  ' 


(13) 


This  new  formulation,  i.e.  Eqs.  (5)-(8),  is  based  on  a  set  of 
coupled  partial  differential  equations  that  track  the  reduced  elec¬ 
trochemical  potential  A*.  reduced  chemical  potential  /x*,  electrical 
potential  0,  current  density  J2,  displacement  current  density  Jdis, 
and  total  current  density JT.  Importantly,  all  the  potential  terms,  A*, 
/i*  and  0,  contain  units  of  voltage.  This  voltage-current  (density) 
formulation  is  especially  convenient  for  electrochemical  modeling 
because  voltage  and  current  signals  can  be  applied  and  measured 
for  the  experiments.  It  also  allows  the  mapping  of  these  differen¬ 
tial  equations  to  an  equivalent  circuit;  thus,  the  physical  meanings 
of  the  charge  transport  can  be  clearly  understood  [6-13].  Finally,  a 
mixed  formulation  from  Eqs.  (l)-(3)  and  (5)-(8),  with  concentra¬ 
tion  Cj,  reduced  electrochemical  potential  A2*>  electrical  potential 
0  and  current  densities  J2,  lead  to  equations  that  are  similar  to 
Newman’s  model  [13,32]. 

For  a  two-charge-carrier  (positive  +  and  negative  -)  particle 
under  electroneutrality  condition,  the  generalized  PNP  equations 
(5)-(8)  under  planar  (m  =  0),  cylindrical  (m  =  l),  and  spherical 
(m  =  2)  symmetry  become  [11,12]  the  following: 


their  coupling.  The  differential  resistor  ( dR+  and  dR-)  and  capacitor 
(dC±em)  elements  are  shown  in  Fig.  1. 

The  surface  voltage  A-(ro)  determines  how  lithium  ions  will  be 
inserted  or  extracted  from  the  electrolyte,  and  the  surface  voltage 
/x*  (L)  determines  how  electrons  will  be  inserted  or  extracted  from 
the  current  collector.  Thus,  the  variables  of  interest  are  the  voltage 
difference  A*  (ro)  -  A+(ro)  and  current  density  J+(r0).  In  addition, 
Eqs.  (14)  and  (15)  can  be  written  as  follows: 

J+-  ^J-  =  ~  £+)•  (18) 

From  the  definition  of  Eq.  (11),  the  difference  of  the  reduced 
electrochemical  potential  is  the  same  as  that  of  the  reduced  chem¬ 
ical  potential 


A*  -  p,*+ =  fl*_  -  fl*+. 


(19) 


These  equations  can  be  simplified  if  two  assumptions  regarding  the 
thermodynamics  and  kinetics  are  made.  First,  it  is  assumed  that  the 
chemical  potential  of  electrons  does  not  depend  on  their  concen¬ 
tration.  Essentially  this  leads  to  the  position  independent  chemical 
potential,  9/x*  /dr  «  0.  Second,  it  is  assumed  that  electronic  mobil¬ 
ity  is  much  larger  than  ionic  mobility;  thus,  electronic  conductivity 
dominates  over  ionic  conductivity  under  the  electroneutrality  con¬ 
dition,  (7+1(7—  ^  0.  With  these  two  assumptions,  Eqs.  (18)  and  (19) 
become  the  following: 


g+  9^+ 
e  dr 


(20) 


,  a  9A+ 

J+  +  dr  ’ 


J-  =  -or- 


9A- 

1F’ 


(14) 

(15) 


Eq.  (20)  can  also  be  written  as  the  following: 

T ms  _  jj  ^0+  —  _  ( 7 +  dc+ 

J+  -  +  dr  -  c^em  dr  ’ 


(21) 


d  (  mi  \  rmrchemd(M-  M+) 

Tr{r  J+)  =  r  c±  at 


Here,  the  definitions  of  Eqs.  (9),  (12)  and  (13)  are  used.  Eqs.  (16) 
and  (17)  lead  to  the  following  equation: 


The  combined  chemical  capacitance  of  the  ion  and  electron  cfiem 
is  the  following: 

=  =  r  *_  n 

Qchem  Qchem  Qchem  e  ' 

These  partial  differential  equations  can  be  mapped  to  an  equiv¬ 
alent  circuit,  as  shown  in  Fig.  1.  Here  the  two-charge  system  is 
a  solid  electrode  particle  with  lithium  ions  (+)  and  electrons  (-). 
The  positions  0  and  r0  represent  the  center  and  surface  of  the 
particle,  respectively.  The  two  resistor  rails  correspond  to  the  trans¬ 
port  of  two  charge  carriers,  and  the  chemical  capacitors  represent 


d_ 

dr 


(rm;^s). 


(22) 


Eqs.  (21)  and  (22)  have  the  same  form  as  the  well-known  Fields 
laws  and  can  be  called  generalized  diffusion  equations.  How¬ 
ever,  the  present  formulations  presents  the  conditions  to  validate 
these  equations.  Second,  it  displays  the  relevant  voltage  variables. 
Third,  it  shows  that  the  chemical  diffusivity  D+  is  related  to  the 
kinetic  (conductivity  cr+)  and  thermodynamic  (chemical  capaci¬ 
tance  C+hem)  properties  and  could  depend  on  the  concentration 
c+. 
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Fig.  2.  (a)  Dimensionless  voltage  -ln[X/(l  -X)]  -g(X-0.5)  for  three  different  thermodynamic  parameters g  =  0,g  =  - 3,  andg=-  5.  B1  and  B2  are  binodal  points  and  SI  and 
S2  are  spinodal  points,  (b)  Dimensionless  chemical  diffusivity  1  +gX(l  -X)  for  three  different  thermodynamic  parameters. 


2.2.  Regular  solution  model 


A  regular  solution  is  a  thermodynamic  model  that  describes  a 
binary  AB  solution  [33,34].  It  is  also  called  the  Frumkin  adsorption 
isotherm  [35-37],  in  which  the  gas  species  and  empty  solid  surface 
sites  can  be  considered  as  A  and  B,  respectively.  In  lithium  intercala¬ 
tion  materials,  lithium  ions  in  the  polyhedral  sites  will  be  denoted 
A  and  empty  intercalation  sites  will  be  denoted  B.  In  the  present 
notation,  the  chemical  potential  of  lithium  ions  is  as  follows: 

M+  =  fil  +  kBTla2^  +kBTg(X-  0.5),  X=^±.  (23) 


The  first  term  on  the  right  hand  side  /x°  is  the  standard  chemical 
potential.  The  second  term  is  the  entropic  contribution,  where  kB 
is  the  Boltzmann  constant  and  T  is  the  absolute  temperature.  The 
third  term  is  the  interaction  term  that  accounts  for  the  interaction 
between  A-A,  B-B  and  A-B  couples,  where  g  is  a  dimensionless 
interaction  parameter.  Variable  c°  is  the  maximum  concentration 
and  X  is  the  dimensionless  concentration.  For  this  regular  solution 
model,  the  volumetric  chemical  capacitance  becomes  the  follow¬ 
ing: 


r'chem 


e2c9  l  l 

- — - 1 - b£ 

kBT  [X  1  -X  5 


(24) 


For  the  kinetics  of  jumps,  the  lithium  ion  mobility  depends  on 
whether  the  potential  jump  site  is  empty,  as  follows: 


b+=b°(  l-X),  (25) 

where  is  a  constant.  This  is  different  from  the  mobility  expres¬ 
sion  of  reference  [38,39],  in  which  the  mobility  is  considered  a 
constant.  The  conductivity  then  becomes  the  following: 


<r+  =  e2c+b+  = 


D°c°e2 

kBT 


X(l-X), 


(26) 


where 


D°  =  kBTb°.  (27) 

Thus,  the  chemical  diffusivity  takes  the  form  of  following  equation: 

D+  =  7^=D°P+SXV-XH-  (28) 

L  + 

As  reported  in  Ref.  [11],  it  is  informative  to  plot  the  dimen¬ 
sionless  voltage  -ln[X/(l -X)] -g(X- 0.5)  in  Eq.  (23)  and  the 


dimensionless  chemical  diffusivity  1  +gX(l  -X)  in  Eq.  (28)  for  dif¬ 
ferent  values  of  interaction  parameters  g.  Three  g  values,  0,  -3,  and 
-5,  are  plotted  in  Fig.  2.  The  Langmuir  adsorption  isotherm  corre¬ 
sponds  to  the  case  ofg= 0.  Whengis  less  than  -4,  both  the  chemical 
capacitance  C^em  and  the  chemical  diffusivity  D+  take  a  negative 
value  at  intermediate  X  values  as  can  be  seen  in  Fig.  2(b).  Because 
chemical  capacitance  is  essentially  the  second  order  derivative  of 
free  energy  with  respect  to  X,  this  system  is  unstable  and  the  initial 
homogeneous  system  will  separate  into  two  phases  [34,40]  when 
g  is  less  than  -4.  In  the  two-phase  coexistence  region,  the  chemi¬ 
cal  potential  is  fixed  so  that  the  equilibrium  thermodynamic  curve 
becomes  the  solid  line  for  g=-5  in  Fig.  2(a).  The  dotted  line  is 
described  by  Eq.  (23 ).  The  two  ends  of  this  dotted  line,  B1  and  B2,  are 
called  binodal  points  and  have  values  of  approximately  0.145  and 
0.855,  respectively.  The  local  minimum  and  maximum  points,  SI 
and  S2,  are  called  spinodal  points  and  are  approximately  0.276  and 
0.724,  respectively.  The  regions  between  B1  and  SI  and  between  S2 
and  B2,  are  metastable.  The  region  between  SI  and  S2  is  the  unsta¬ 
ble  spinodal  region.  Wheng  is  equal  to  -4,  the  chemical  capacitance 
goes  to  infinity,  and  the  chemical  diffusivity  goes  to  0  at  X=0.5 
(Fig.  2(b)),  which  was  discussed  in  Ref.  [11]. 

In  this  study,  the  regular  solution  model  with  g= 0,  g=  -  3,  and 
g=-  5  was  chosen  to  represent  a  solid  solution  with  constant 
diffusivity,  a  solid  solution  with  variable  diffusivity,  and  a  phase 
transformation  material,  respectively.  During  phase  transforma¬ 
tion,  a  phase  boundary  develops  and  moves  between  coexisting 
phases.  In  the  present  work,  a  diffuse  interface  model  using  the 
phase  field  method  was  utilized  to  study  the  movement  of  the 
phase  boundary.  In  addition  to  its  widespread  application  in  solidi¬ 
fication  and  morphology  evolution  [41  -44],  the  phase  field  method 
has  recently  been  applied  in  a  few  works  to  the  study  of  electro¬ 
chemistry  in  the  context  of  electrodeposition  [45,46]  and  phase 
transformation  in  battery  electrodes  [38,39,47-49]. 

In  the  phase  field  method,  the  chemical  potential  [39]  is  repre¬ 
sented  by  the  following: 

M+  =  M+  +  kBT In  2L^  +  kBTg(X  -  0.5) 

(29) 

Compared  to  Eq.  (23),  the  extra  fourth  term  is  the  Cahn-FIilliard 
interfacial  energy  term  [50],  where  y  is  the  dimensionless  inter¬ 
facial  energy  parameter.  This  term  goes  to  0  in  the  case  of  solid 
solutions  with  g  =  0  and  g=  -  3,  i.e.  Eq.  (23).  The  term  y  is  taken  to 
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Fig.  3.  Schematic  depiction  of  interfacial  ionic  transport  at  the  solid  |  liquid  interface. 
The  relevant  concentrations  c+  and  reduced  electrochemical  potentials  jl*+  of  the 
solid  and  liquid  are  labeled  along  with  the  interfacial  length  Ssl- 


be  10-5  in  this  work.  This  interfacial  parameter  is  related  to  the 
interface  between  two  coexisting  phases  in  the  solid,  but  not  to  the 
solid|liquid  interface  to  be  discussed  shortly. 


2.3.  Boundary  conditions 


If  the  solid  electrode  particle  is  in  contact  with  the  liquid  elec¬ 
trolyte,  the  transport  of  lithium  ions  occurs  across  the  solid  [liquid 
interface.  The  latter  is  shown  schematically  in  Fig.  3, where  the 
interfacial  length  is  ^.Variables  c+(r0)  and  c+(L)  are  the  lithium 
ion  concentrations  in  the  solid  and  liquid,  respectively,  on  the  two 
sides  of  the  interface.  The  driving  force  is  the  difference  between 
the  reduced  electrochemical  potentials  of  lithium  ions,  A+(ro)  and 
£*  (I).  The  jump  is  assumed  to  exponentially  depend  on  the  driv¬ 
ing  force  as  an  Arrhenius  rate  expression.  The  jumping  probability 
also  depends  on  the  direction.  Specifically,  jumping  from  liquid  to 
solid,  the  jump  rate  is  proportional  to  the  fraction  of  available  sites 
given  by  1  -X(r0),  while  the  jump  from  solid  to  liquid  has  no  such 
restriction. 

Thus,  the  current  density  from  solid  to  liquid  and  liquid  to  solid 
can  be  written  as  follows: 


Jsl  =  X(r0)exp 

°SL 


£+(r0) -£+(£) 


kBT/e 


Jls  =  -*(r0)]exp 

OSL  CV 


£+(i)-£+(r0) 

kBT/e 


(30) 

(31) 


where  ks  and  kL  are  parameters  that  incorporate  the  jumping  prob¬ 
abilities.  Combining  bidirectional  fluxes,  the  current  density  across 
this  interface  from  solid  to  liquid  can  be  written  as  follows: 


J+(r o)  =  ^-X(r0)exp 

OSL 


£+(ro)  -  £+(!•) 


kBT/e 


^^}[i_X(r0)]exp 


8sl 


£+(r0) -£+(£) 

kBT/e 


(32) 


Clearly,  the  application  of  the  above  rationale  to  the  jump  from 
position  x  to  x  +  dx  in  the  solid  provides  the  following  equation: 


J+  =  ^X(x)[l  -X(x  +  dx)]  exp 


£*(x)-  jL\{x  +  dx) 


k(x  +  dx) 
dx 


X(x  +  dx)[l  -  X(x)]  exp 


kBT/e 
£+(x)  -  £+(x  +  dx) 

W~e 


(33) 


As  dx  goes  to  0,  Eq.  (33)  can  be  written  as  the  following: 


J+  =  — 2fc(r)X(r)[l -X(r)]^jP-.  (34) 

This  expression  is  consistent  with  Eqs.  (5)  and  (26),  as  shown 
before. 

At  the  solid|liquid  interface,  the  expression  of  Eq.  (32)  can  be 
simplified  in  different  ways.  First,  it  can  be  written  as  the  following: 


J+(ro)  = 


=  Jo| 


exp 


£+(r0)-£^(L) 


kBT/e 


-  exp 


-a 


£+(ro) -£+(£) 


kBT/e 


=  2J0  sinh 


£+(ro) -£+(!•) 
“  kBT/e 


(35) 


This  equation  has  a  form  similar  to  that  of  the  Butler-Volmer  equa¬ 
tion  [14,51-53].  The  apparent  exchange  current  density  J0  and 
symmetry  factor  a  can  be  represented  as  the  following: 


q.  1  We  ln  ksX(r0) 

2  £*(r0)-£*(L)  kL[\  -X(r0)]c+(L)/c°' 


(36) 

(37) 


This  formulation  is  also  similar  to  the  symmetrized  form  of  the 
transport  equation  under  a  large  gradient  driving  force  [54]. 
Because  the  classical  Butler-Volmer  equation  can  be  approximated 
by  the  Tafel  and  linearized  equations  [14],  two  similarly  simpli¬ 
fied  forms  of  the  generalized  Butler-Volmer  equation,  Eq.  (35),  can 
be  obtained  as  de  Donder’s  equation  and  the  resistive  equation, 
respectively.  First,  if  the  argument  in  the  hyperbolic  sine  function 
is  large  and  positive,  Eq.  (35)  can  be  written  as  follows: 

J+(ro)  =Jo  jexp 

This  equation  has  the  form  of  de  Donder’s  equation,  as  in  Ref. 
[47].  The  term  of  1  in  the  curly  bracket  can  be  ignored  because  it  is 
much  smaller  than  the  exponential  term  and  it  provides  the  Tafel 
form  [14]. 

Second,  if  the  term  in  the  hyperbolic  sine  function  of  Eq.  (35)  is 
small,  it  can  be  linearized  as  follows: 


£+(fo) -£+(£) 

?  kBT/e 


- 1 


(38) 


J+(fo)  = 


2aJ0 

kBT/e 


[£;(r0) -£;a)]. 


(39) 


This  equation  can  be  called  the  generalized  Chang-Jaffe  equation. 
In  the  conventional  Chang-Jaffe  equation  [55],  the  term  in  front  of 
the  square  bracket  is  a  constant,  and  the  difference  term  inside  the 
square  bracket  is  the  concentration  difference.  Furthermore,  if  the 
term  before  the  square  bracket  is  assumed  to  be  a  constant,  Eq.  (39) 
becomes  the  following: 


(40) 


where  R-mt  is  a  constant.  This  expression  was  used  in  the  study 
of  interfaces  involving  mixed  conductors  [10,56,57]  and  can  be 
called  a  resistive/ohmic  boundary  condition.  This  was  also  called 
the  cell-impedance-controlled  boundary  condition  [58].  When  Rint 
is  small,  it  is  expected  that  variable  /z+(ro) -  £+(0  will  also  be 
small;  this  can  be  called  the  reversible  boundary  condition.  In  this 
case,  J+(r0)  can  be  determined  by  the  “diffusion  equations”,  i.e. 
Eqs.  (20)  and  (22),  which  relate  J+  to  £+(r0);  this  can  be  called 
a  diffusion-controlled  mechanism.  When  Rint  is  large,  it  can  be 
called  an  interface-controlled  mechanism  or  kinetics-controlled 
with  ohmic  behavior.  For  intermediate  values  of  Rint,  it  can  be  called 
a  mixed-controlled  mechanism  [26]. 

Eq.  (40)  is  the  interfacial  solid  [liquid  boundary  condition  that  is 
employed  in  the  present  study  and  is  equivalent  to  replacing  the 
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electrolyte  box  with  a  resistor  in  Fig.  1.  At  the  other  boundary,  i.e. 
the  particle  center,  the  flux  is  zero  due  to  symmetry,  as  is  described 
below: 


J+(  0)  =  0.  (41) 

Finally,  in  Fig.  1,  the  voltage  difference  /z*  (r0)  -  was  used 
as  the  external  voltage  Vext  and  the  current  density  J+(r0)  was  used 
as  the  external  current.  Thus,  Eq.  (40)  was  transformed  to  become 
the  following: 


Then,  Eq.  (45)  becomes  the  following: 


‘"If  -  s<-‘"« 


'"T^x+*(x-05)-4s 


m = o 

i/r(\,r  =  0)-f{\)  +  Av 
J(  )  R 


(48) 


,  ,  .  H*+{ro )  -  H-(r0)  +  Vext 
J+{r0)  = - „ - 


(42) 


The  equivalence  of  reduced  electrochemical  potential  difference 
to  reduced  chemical  potential  difference,  Eq.  (19),  was  once  again 
applied.  In  addition,  the  initial  conditions  were  always  set  to  the 
open  circuit,  where  the  current  was  0  and  the  voltage  was  the  open 
circuit  voltage  Vext(0): 


WO)  =  (r0,  t  =  0)  -  /4(r0,  t  =  0)  (43) 

The  voltage  difference  AV=Vext  -  Vext{0)  after  the  current  or  volt¬ 
age  perturbation  obeyed  the  following  equation: 


,  ,  .  M+(»o)  -  M+fo.  t  =  0)  +  AV 

J+lro)  =  - 5 - 


(44) 


Flere,  the  chemical  potential  of  electrons  is  a  constant  as  mentioned 
before. 

For  convenience,  the  partial  differential  equations,  Eqs.  (20)  and 
(22),  along  with  the  boundary  conditions,  Eqs.  (41)  and  (44),  and 
expressions  of  relevant  parameters,  Eqs.  (26)  and  (29),  are  listed 
below: 


a+  3^; 


jms  _  _  ~-t~  ■ 

J+  “  e  dr 
rm°£=-JL(rmjms) 


3 1 

J™(0)  =  0 


JT(r o)  = 


M+(fo)  ~  M+(rp,  t  =  0)  +  AV 

^int 


= 


£)0  C°  6 2 

U+C+e-X(\-X) 


kBT 


X 


=  +  kBT  In  y— —  +  kBTg(X  -  0.5) 


1  d 


-^VaFraF 


,ax 


(45) 


2.4.  Dimensionless  form 


where  ^(1,  t  =  0)  +  Av  is  the  voltage  in  the  liquid  and  ^(1)  is  the 
voltage  on  the  solid  surface.  The  goal  of  the  present  work  is  to  find 
the  relation  between  j(l)  and  ^(1,  r  =  0)  +  Av. 

In  the  battery  field,  a  commonly  used  notation  for  current  is  the  C 
rate.  C/n  is  the  current  that  takes  n  hours  to  fully  charge  or  discharge 
a  battery.  For  a  single  particle  with  maximum  concentration  c+, 
the  mass  flux  density  corresponding  to  the  C/n  rate  becomes  the 
following: 


THIS 

'+ 


(m  +  1  )n  •  3600 ' 


(49) 


Then,  the  dimensionless  mass  flux  density  in  C  rate  is  described  as 
the  following: 


/=— 1 _ *5_. 

1  m  +  1  3600n 


(50) 


For  example,  a  dimensionless  flux  of  1  through  the  surface  of  a 
planar  (m  =  0)  particle  corresponds  to  the  C  or  C/1 0  rates,  when  the 
characteristic  time  to  is  1  h  or  10  h,  respectively. 


2.4.1.  Phase  transformation  material 

For  the  phase  transformation  material,  Eq.  (48)  becomes  the 
following: 


f  = 


lnT^x+s(x_a5) 


1  3 


Yxm  dx 


(51) 


,3X 

9r 


d_ 

dx 


-XmX(l  -X) 


dxf/ 

dx 


(52) 


Eqs.  (51 )  and  (52)  are  two  coupled  second-order  partial  differential 
equations  on X  and  \J/.  For  Eq.  (51 ),  the  boundary  conditions  are  the 
following: 

ax(0)  3X(1)  _ 

y^r  =  y^r=°-  (53) 

This  indicates  that  phase  transformation  will  not  occur  (zero  inter¬ 
facial  energy  between  two  phases)  at  the  center  and  at  the  surface 
of  the  particle  [47].  The  boundary  conditions  for  Eq.  (52)  are  the 
last  two  equations  in  Eq.  (48). 


For  numerical  simulations,  it  is  convenient  to  use  the  following 
dimensionless  variables: 


r 


t 

£o  ’ 


j  = 


jms 


C%/t0  ’ 


ty  = 


M+-M+ 

kBT 


2.4.2.  Solid  solution  material 

For  the  solid  solution  material,  the  second  order  partial  differ¬ 
ential  equation  Eq.  (51)  becomes  the  following  expression: 


\[r  =  — 


In 


X 

1  -X 


+  g(X  -  0.5) 


(54) 


R  = 


kBTr0/(D°c°e2)’ 


Av  = 


eAV 


(46) 


Because  it  does  not  involve  partial  derivatives,  Eq.  (52)  becomes  the 
following: 


The  characteristic  time  t0  is  represented  by  the  following: 


TdX_d_ 
dr  dx 


dx 

xm[l+gX(l-X)]g 


(55) 


(47)  This  is  the  same  second  order  partial  differential  equation  used  in 
the  present  author’s  previous  work  [11]. 
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Table  1 

Partial  differential  equations  along  with  parameters  used  in  thermodynamics,  initial  conditions  (IC),  and  boundary  conditions  (BC). 


II 

>< 

1 

>< 

sjj 

vm  3X  3  r  vm;) 

*  3?  -  3^-X  J) 

*  =  - 

Solid  solution 
g= 0,  Y  =  0 

Solid  solution 
g=~ 3,  Y  =  0 

Phase  transformation 
g=  —  5,  y  =  10“5 

Step  current 

IC 

X0  =  0.99 

u= 

=  0)  =  -ln[Xo/(l-X0)]-g(X0-0.5) 

BC 

V3X(0)  _  V3X(1)  _  n 

y  dx  ~  y  dx  — 

j(0)  =  0 

1(1)  =  0.05,  0.1,  0.2,  0.5 

Step  voltage 

IC 

Xo  =  0.9 

=  0)  =  -ln[Xo/(l-X0)]-g(X0-0.5) 

BC 

v  3X(0)  _  3X(1)  _  n 

y  dx  —  y  dx  - u 

1(0) =0 

j(X)  —  Vf(i>r=o)-V''(i)+Av 

Av  —  0.5,  1,  2,  5]R  =  0A,  1,  10 

Linear  sweep  voltage 

IC 

X0  =  0.99 

yr(l,  r  = 

=  0)  =  -  ln[X0/(l  -  Xo)]  -g(X o  -  0.5) 

BC 

v  3X(0)  _  V3X(1)  _  n 

y  dx  ~  y  dx  — 

1(0) =0 

—  0.05,  0.1,  0.2,  0.5,  1; 

R  =  0A,  1,  10 

Sinusoidal  current 

IC 

Xo  =  0.9 

l/r(l,  T 

II 

2 

II 

l 

3 

1 

(X? 

£ 

1 

O 

Ln 

BC 

v  3X(0)  _  V3X(1)  _  n 

y  dx  -  y  dx  - u 

1(0) =0 

X 1 )  =lacsin(27r/r ), 

jAC= o.oooi,  o.ooi,  o.oi,  0.05; 

/=0.05-1000 

Unless  otherwise  specified,  the  initial  condition  is  chosen  as  a 
constant  concentration  profile  as  follows: 

X(x,  r  =  0)=X0,  (56) 


A  small  perturbation  of  any  form  of  j(l)  will  likely  cause  small 
changes  in  X  and  x/s.  If  both  dx/z/dX  and  X(1  -X)  are  assumed  to  be 
constant  during  the  course  of  perturbation  and  initial  concentration 
profile  is  uniform X0,  perturbations  of  Eq.  (58)  yield  the  following: 


x/z(x,  r  =  0)  = 


lnT^+g(Xo“°'5) 


(57) 


The  above  partial  differential  equations  are  solved  with  Comsol 
Multiphysics  software  [59].  The  dimensionless  spatial  grid  is  uni¬ 
form,  fixed  in  time,  and  composed  by  500  segments.  It  was  found 
the  solution  converges  between  500  and  5000  elements,  with  the 
exception  of  some  cases  in  impedance  spectroscopy  (details  to  be 
discussed  later).  The  equations  are  solved  using  the  finite  element 
method  with  a  Lagrange  quadratic  basis,  ensuring  quadratic  con¬ 
vergence  of  the  spatial  discretization.  In  time  the  problem  is  solved 
using  implicit  finite  backward  finite  differences  of  variable  order. 


3.  Results  and  discussion 


Before  presenting  the  numerical  results  of  Eq.  (48),  the  known 
small-signal  analytical  solutions  [  1 2,30]  are  discussed.  The  first  two 
partial  differential  equations  in  Eq.  (48)  can  be  formally  rewritten 
as  the  following: 


dxdi 

dx//  dr 


d_ 

dx 


-xmX(l  -X) 


dx/f 

dx 


(58) 


,  dAxfs 
dx 


dxfr 

ax 


Xo(l-Xo) 

xo 


=  D(X o) 


d_ 

dx 


(59) 


where  the  change  of  voltage  A  ^  is  the  variable  of  interest  and  D(X0 ) 
has  the  following  form: 

D(Xo)  =  l+gX0(l-X0).  (60) 


After  the  Laplace  transformation  with  boundary  condition  j(0)  =  0 
[12,30],  in  the  planar  case  (m  =  0),  Eq.  (59)  can  be  solved  to  provide 
the  following: 

Ay<l)(s)  1  cosh  yf  s/D(X0) 

ms)  Xo(l-Xo)  -^/s/D(X0) 

where  the  overbar  represents  the  Laplace  transformed  variable, 
and  s  is  a  complex  argument.  The  time-domain  current/voltage 
perturbation  signal,  e.g.  current  step  or  voltage  step,  can  be  trans¬ 
formed  to  the  Laplace  domain,  and  then,  Eq.  (61)  can  be  inversely 
transformed  back  to  the  time  domain  to  obtain  the  corresponding 


9=0 
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Fig.  4.  Evolution  of  the  concentration  profile  under  a  constant  charging  current  j(l )  of  0.05  for  three  different  thermodynamic  parameters  g= 0,  g=  -  3,  and  g=  -  5.  Dimen¬ 
sionless  times  are  shown  on  each  curve.  Bl,  B2,  SI,  and  S2  points  are  the  binodal  and  spinodal  points  discussed  in  Section  2.2.  Particle  surface  and  center  are  at  the  position 
of  1  and  0,  respectively. 
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voltage/current  signal.  When  s  =  ioo ,  where  i  is  the  imaginary  num¬ 
ber  and  oo  is  the  angular  frequency,  the  Laplace  transform  becomes 
the  Fourier  transform  and  Eq.  (61)  provides  the  impedance.  Short- 
time  and  long-time  analytical  expressions  obtained  from  Eq.  (61) 
for  the  current  step  (Section  3.1),  voltage  step  (Section  3.2), 
and  sinusoidal  current  (Section  3.4)  signals  were  compared  with 
numerical  results  for  signals  of  an  arbitrary  magnitude. 

For  the  ease  of  comparison,  parameters  of  thermodynamic 
parameters,  initial  conditions  (IC),  and  boundary  conditions  (BC) 
used  in  the  various  techniques  of  Sections  3.1 -3.4,  are  listed  in 
Table  1,  along  with  partial  different  equations  that  were  solved. 


region,  e.g.  at  r  =  3,  before  it  reached  the  unstable  region.  As  can 
be  seen  from  Fig.  2(a),  the  concentration  in  the  metastable  region 
corresponded  to  a  voltage  overshoot  beyond  the  voltage  plateau. 
This  overshoot  has  been  experimentally  observed  for  LiCo02  as 
it  goes  through  the  insulator-metal  or  03-1  to  03-11  phase  tran¬ 
sition  [60-62]  and  for  the  triphylite-heterosite  phase  transition  of 
LiFeP04  [63]. 

Finally,  the  small-signal  solution  under  the  context  of  GITT  are 
discussed.  For  the  constant  current  j(l),  the  short-time  voltage 
response  can  be  obtained  from  the  Laplace-inverse  Laplace  trans¬ 
formation  of  Eq.  (61)  [12,30]: 


3. t.  Constant- current  chronopotentiometry 

In  the  constant-current  chronopotentiometry  study,  a  constant 
current  j(l)  was  applied  and  the  voltage  was  monitored  as  a  func¬ 
tion  of  time.  As  discussed  in  Section  2.4,  for  the  present  single 
particle,  voltage  in  the  liquid  ^(1,  r  =  0)  +  Av  was  the  variable 
of  interest.  Flowever,  because  the  solid  surface  was  connected  to 
the  liquid  by  a  resistor,  the  extra  voltage  difference  could  be  easily 
obtained  by  multiplying  resistance  with  the  constant  current.  Thus, 
the  voltage  at  the  particle  surface  i/s(  1 )  was  studied  instead  and  this 
also  applied  to  the  controlled-current  impedance  spectroscopy  in 
Section  3.4. 

Initially,  the  concentration  X0  was  a  uniform  value  of  0.99  in  the 
solid,  and  the  simulation  was  performed  in  the  charging  process 
for  a  cathode  particle.  The  current  values  j(l)  were  0.05,  0.1,  0.2, 
and  0.5  in  dimensionless  form.  The  time  stepping  ended  when  the 
surface  concentration  X(l)  reached  0.01.  The  time  step  was  0.005 
for  a  current  value  of  0.5;  otherwise,  the  time  step  was  0.02. 

The  evolution  of  the  concentration  profiles  in  the  solid  under 
a  constant  current  of  0.05  is  shown  in  Fig.  4.  In  the  case  of  solid 
solutions  with  g=0,  -3,  the  concentration  on  the  surface  x=l 
decreased  under  the  application  of  a  constant  charging  current  and 
gradually  diffused  into  the  center.  In  the  case  of  the  phase  trans¬ 
formation  material  with  g=  -  5,  the  area  defined  by  the  two  inner 
horizontal  dotted  lines  (SI  and  S2)  represents  the  unstable  spin- 
odal  region,  as  discussed  in  Section  2.2.  The  two  shells  formed  by 
the  inner  dotted  lines  and  outer  dotted  lines  (B1  and  SI ;  S2  and  B2) 
are  the  metastable  regions.  Clearly,  the  evolution  of  the  concentra¬ 
tion  profile  was  similar  to  that  in  the  solid  solution  before  reaching 
the  unstable  region  X  =  0.724  ( r  <  4).  Then,  phase  separation  started 
and  the  phase  boundary  moved  towards  the  center  of  the  particle 
(4  <  r  <  16).  Finally,  the  particle  was  in  another  single-phase  region 
for  r  >  18.5. 

It  was  interesting  to  study  how  fast  the  phase  boundary  moved 
and  how  this  process  depended  on  the  current.  The  position  and 
velocity  of  the  phase  boundary  at  different  currents  (0.05,  0.1,  0.2, 
and  0.5)  are  shown  in  Fig.  5.  First,  Fig.  5(a)  indicates  that  it  moved 
linearly  with  time  for  the  current  values  investigated: 

XPB  +  vPBr>  (62) 

where  x£g  and  v ££  are  the  phase  boundary  position  and  veloc¬ 
ity  respectively.  Variable  XqC  is  a  constant.  Fitting  the  x£g  data  to 
Eq.  (62)  yielded  v ££,  and  it  was  plotted  as  a  function  of  current  in 
Fig.  5(b).  The  phase  boundary  moved  faster  under  higher  current, 
as  intuitively  expected. 

The  capacity-voltage  results  are  shown  in  Fig.  6.  The  capacity, 
the  product  of  current  j(l)  and  time  r,  was  used  to  fit  the  battery 
field  convention.  The  thermodynamic  curves  in  Fig.  2  are  also  plot¬ 
ted  in  dashed  lines  as  a  reference.  As  the  applied  current  increased 
from  0,  polarization  increased  for  all  three  thermodynamic  cases, 
as  expected.  In  addition,  the  voltage  overshoot  was  observed  at  the 
beginning  of  charging  process  in  the  phase  transformation  case. 
Fig.  4  shows  that  surface  concentration  entered  the  metastable 


AtA(1)  _  2/(1)  /d(X o) 

Sc  X0(1-X0)V 


(63) 


where  the  right-hand  side  is  a  constant.  This  is  similar  to  the  expres¬ 
sion  in  the  GITT  Ref.  [17].  Results  of  plotting  short-time  (r  <  1 )  data 
of  Fig.  6  in  the  form  of  Eq.  (63)  are  shown  in  Fig.  7.  This  equation 
did  not  hold  even  for  the  lowest  current  of  0.05,  which  suggested 
that  the  small-signal  approximation  is  not  satisfied  at  this  current. 
However,  this  approximation  is  valid  for  a  different  initial  condition 
that  will  be  discussed  in  Section  3.5. 


3.2.  Potential  step  chronoamperometry 

In  the  potential  step  chronoamperometry,  a  constant  voltage 
step  Av  was  applied  and  the  current  j(l )  was  monitored  as  a  func¬ 
tion  of  time.  The  boundary  condition  for  j(l)  in  Eq.  (48)  and  three 
R  values,  0.1,  1,  and  10,  were  studied.  These  values  were  chosen 
to  study  different  mechanisms,  such  as  diffusion  control,  mixed 
control,  and  interface  control. 

Initially,  the  concentration  X0  was  a  uniform  value  of  0.9  in  the 
solid,  and  the  simulation  was  performed  in  the  charging  process. 
The  voltage  steps  Av  were  0.2,  0.5, 1,  2,  and  5.  The  time  was  incre¬ 
mented  in  steps  0.02  tillj(l)  fell  below  0.001. 

Similar  to  the  case  of  constant-current  chronopotentiometry, 
the  evolution  of  concentration  profiles  was  studied  first.  Fig.  8 
shows  an  example  of  the  case  with  Av  =  0.5  and  R  =  l.  If  there 
was  no  solid |  liquid  interface,  the  voltage  step  implied  that  sur¬ 
face  concentration  was  fixed,  as  in  Ref.  [11].  With  an  interface,  the 
voltage  step  translated  to  a  variable  concentration  step.  Overall, 
the  concentration  steps  led  to  steeper  concentration  profiles  com¬ 
pared  to  the  current  steps,  in  solid  solutions  ofg  =  0,  -3  in  Fig.  4. 
In  the  case  of  phase  transformation  with  g=  -  5,  the  change  in  the 
concentration  profile  under  constant  potential  was  similar  to  the 
case  of  constant  current  in  Fig.  4.  A  voltage  step  of  0.5  eventually 
changed  the  surface  from  the  high  concentration  phase  to  a  low 
concentration  phase,  which  was  accompanied  by  the  phase  bound¬ 
ary  movement  towards  the  center.  For  the  smallest  voltage  step  of 
0.2,  no  phase  transformation  was  initiated. 

Again,  the  positions  of  the  phase  boundary  were  tracked 
(Fig.  9(a))  with  different  voltage  step  Ay  values  and  interfacial  resis¬ 
tance  R  values.  The  phase  boundary  moved  linearly  with  the  square 
root  of  time: 

xpb  ^xoV  +  i'pfiT0-5,  (64) 

where  x^  and  v ^  are  phase  boundary  position  and  velocity  con¬ 
stant,  respectively.  Variable  x^v  is  a  constant.  As  the  interfacial 
resistance  was  increased,  e.g.  R  =  10,  the  square  root  dependence 
became  less  obvious.  Fitting  the  x£g  data  to  Eq.  (64)  yielded  , 
and  it  was  plotted  as  a  function  of  the  voltage  step  in  Fig.  9(b).  The 
phase  boundary  moved  faster  under  larger  voltage  steps,  as  intu¬ 
itively  expected.  In  addition,  the  phase  boundary  velocity  showed 
little  variation  between  R  =  0.1  and  R  =  1,  which  suggested  that 
both  cases  were  diffusion-controlled.  It  can  be  expected  that  R  =  10 
would  correspond  to  a  mixed-controlled  mechanism,  and  larger 
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Fig.  5.  (a)  Phase  boundary  positions,  and  (b)  phase  boundary  velocities  under  different  values  of  constant  current  j(l)  (0.05,  0.1,  0.2,  and  0.5)  for  the  phase  transformation 
material  with  g=  -  5. 


g=0  g=-3  g=-5 


Capacity 

Fig.  6.  Voltage  \J/(  1)  as  a  function  of  capacity  under  three  different  currents  j(l)  (0.05,  0.2  and  0.5)  for  three  different  thermodynamic  parameters  g=0,  g=  -  3,  and  g=  -  5. 
Dashed  lines  represent  the  thermodynamic  voltage  curve  in  Fig.  2(a).  The  voltage  overshoot  can  be  clearly  observed  in  the  inset  of  g=  -  5. 


R  values  would  lead  to  an  interface-controlled  mechanism.  As  R 
became  larger,  e.g.  100,  the  linear  time  dependence  was  observed 
at  long  time  periods  (not  shown).  In  the  sharp-interface  moving 
boundary  model,  the  square  root  and  linear  dependences  were  also 
obtained  for  diffusion  control  and  interface  control  mechanisms, 
respectively  [64]. 


Finally,  the  relation  between  the  voltage  step  Av  and  current 
response  j(l)  are  shown.  First,  small-signal  analytical  expressions 
from  Eq.  (61)  are  discussed.  If  no  interface  effect  was  present,  Av 
was  equivalent  to  A^(l).  The  short-time  (Eq.  (65))  and  long-time 
(Eq.  (66))  responses  of  the  current  j(l)  were  obtained  from  the 


g=0  g=-3  g=-5 


Fig.  7.  Short-time  characteristic  curves  according  to  Eq.  (63)  under  three  different  currents  j(l )  (0.05, 0.2  and  0.5)  for  three  different  thermodynamic  parameters  g= 0,  g=  -  3, 
and  g  =  -  5.  Dashed  lines  are  the  plateaus  expected  from  the  small-signal  assumption. 
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Fig.  8.  Evolution  of  the  concentration  profile  under  a  constant  voltage  step  Av  =  0.5  and  an  interfacial  resistance  R  =  1  for  three  different  thermodynamic  parameters  g  =  0, 
g=-  3,  andg=-  5.  Dimensionless  times  are  shown  on  each  curve.  Bl,  B2,  SI,  and  S2  points  are  the  binodal  and  spinodal  points  discussed  in  Section  2.2.  Particle  surface  and 
center  are  at  the  position  of  1  and  0,  respectively. 


Laplace-inverse  Laplace  transform  [12,30]  as  follows: 

(65) 

\AtD(X0) 


j(l)  =  2X0(1  -X0)exp 


~D(X0)t 


A^r(l), 


(66) 


where  A^(l )  is  the  voltage  step.  A  slightly  different  form  of  Eq.  (66) 
was  also  presented  in  Ref.  [12]  using  a  different  series  expansion. 


Eq.  (65)  suggests  that  the  product  of  current  and  the  square  root 
of  time  is  a  constant,  also  known  as  the  Cottrell  equation  [14].  Eq. 
(66)  suggests  that  the  logarithm  of  current  changes  linearly  with  r 
at  long  time  periods.  These  are  essentially  the  same  equations  as  in 
the  PITT  Ref.  [18]. 

The  numerical  simulation  in  this  work  involved  the  solid  [liquid 
interface  and  made  no  assumption  of  small  signals.  Under  these  cir¬ 
cumstances,  the  analytical  solutions  were  not  available.  However, 
the  results  of  the  solid  solutions  are  still  presented  in  the  form  of 


Fig.  9.  (a)  Phase  boundary  positions,  and  (b)  phase  boundary  velocity  constant  under  different  values  of  voltage  steps  Av  (0.5,  1,  2,  and  5)  and  interfacial  resistances  R  for 
the  phase  transformation  material  withg=  -  5. 
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Eqs.  (65)  and  (66)  for  comparison  as  dashed  lines  in  Fig.  10(a)-(d), 
while  the  regular  current  vs.  time  plot  is  displayed  in  Fig.  10(e)  for 
the  phase  transformation  material. 

For  the  short-time  response  of  the  solid  solutions,  Fig.  10(a)  for 
g  =  0  and  Fig.  10(b)  for  g=-3,  the  Cottrell  behavior  can  only  be 
observed  at  small  R  =  0.1  values  and  small  voltage  step  of  Av  =  0.2. 
A  small  value  of  R  suggested  that  it  was  under  diffusion  control, 
and  a  small  voltage  step  Av  implied  a  small-signal  perturbation. 
These  two  conditions  led  to  results  that  were  similar  to  those  in 
Eqs.  (65)  and  (66).  As  R  became  larger,  i.e.  toward  to  interface 
control,  or  Av  became  larger,  i.e.  large  signal,  peak  shaped  plots 
were  observed  in  Fig.  10(a)  and  (b).  These  types  of  peak  shapes 


have  been  observed  experimentally  for  graphite  [65]  and  LiFeP04 

[66]. 

For  the  long-time  response  of  solid  solutions,  Fig.  10(c)  forg=0 
and  Fig.  10(d)  for  g=-3,  a  linear  behavior  can  be  observed  in  all 
the  R  and  Av  values  investigated.  However,  the  slope  only  cor¬ 
responded  to  values  predicted  from  Eq.  (66)  in  situations  where 
approximations  of  diffusion  control  and  small-signal  can  be  satis¬ 
fied,  such  as  low  R  and  Av  values.  In  the  case  of  g=  -  3,  Fig.  10(d), 
deviations  of  predicted  slopes  can  be  observed  for  most  values 
studied. 

For  the  case  of  phase  transformation  with  g=  -  5,  the  current 
response  following  the  voltage  step,  Fig.  10(e),  exhibited  some 


Time 


Time 


R=10  R=1  R=0.1 


Fig.  10.  Short-time  characteristic  curves  according  to  Eq.  (65)  for  (a)g  =  0  and  (b)g  =  -3  under  three  different  voltage  steps  Av  (0.2,  0.5,  and  1)  and  interfacial  resistance  R 
values.  Dashed  lines  represent  the  expected  plateau  from  small-signal  assumption.  Long-time  characteristic  curves  according  to  Eq.  (66)  for  (c)  g=0  and  (d)  g=-3  under 
three  different  voltage  step  and  interfacial  resistance  values.  Dashed  lines  represent  the  expected  straight  lines  from  small-signal  assumption,  (e)  Current  response  for  phase 
transformation  material  g=  -  5  under  three  different  voltage  step  and  interfacial  resistance  values. 
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Fig.  10.  (  Continued ). 


interesting  behavior.  Initially,  the  material  was  in  the  single-phase 
region.  When  the  voltage  step  Av  was  small,  e.g.  0.2,  the  material 
stayed  as  a  single  phase  for  all  three  R  values,  thus,  the  current 
decayed  in  a  manner  similar  to  Eqs.  (65)  and  (66).  As  Av  increased 
to  0.5,  the  current  spike  was  observed  for  R  values  of  10  and  1.  The 
comparison  of  the  current  response  in  Fig.  10(e)  and  concentra¬ 
tion  evolution  in  Fig.  8  for  the  same  set  of  data  ( Av  =  0.5  and  R  =  1 ) 
suggested  that  the  current  spike  corresponded  to  the  initiation  of 
phase  separation,  and  the  kink  at  approximately  r  =  1 5  marked  the 
finish  of  the  separation.  This  behavior  was  observed  more  clearly 
with  values  of  Av  =  0.5  and  R  =  10,  which  showed  the  initial  single¬ 
phase  decay,  initiation  of  phase  transformation,  two-phase  decay, 
and  finally,  another  single-phase  decay.  The  two  separate  single¬ 
phase  decaying  stages  were  similar  to  those  of  solid  solutions  with 
exponential  decay,  Fig.  10(c)  and  (d).  The  two-phase  decay  corre¬ 
sponded  to  the  phase  boundary  movement,  which  appeared  to  be 
slower.  As  Av  increased  further  to  1  (R  =  1 ),  the  initial  spike  disap¬ 
peared,  but  the  kink  was  still  visible.  This  situation  still  progressed 
through  the  phase  transformation  cycle  as  in  the  case  of  Av  =  0.5. 
Physically,  the  initiation  of  phase  separation  and  the  phase  bound¬ 
ary  movement  corresponded  to  nucleation  and  growth  of  the  new 
phase. 

First,  the  initial  current  increase  upon  voltage  step  was  observed 
experimentally  for  graphite  [67,68]  and  LiFeP04  [66],  though  in  the 
form  of  a  hump  instead  of  a  spike.  The  hump  was  explained  by  the 
nucleation  of  new  phases  [66-68]  and  was  quantitatively  modeled 
by  the  Johnson-Mehl-Avrami  equation  [49,69].  Because  the  cur¬ 
rent  spike  in  Fig.  10(e)  was  associated  with  the  initiation  of  phase 
separation,  it  was  expected  that  a  homogenized  model  of  many  par¬ 
ticles  with  a  distribution  of  particle  sizes  and  shapes  would  lead 
to  a  dispersion  of  time  constants  to  form  the  hump  observed  in 


those  experimental  composite  electrodes.  Second,  the  multi-staged 
current  decay  was  also  obtained  from  the  sharp-interface  moving 
boundary  problem  [58]  and  has  been  experimentally  observed  for 
Li4Ti50i2  and  Li§V205  [70].  Finally,  the  square  root  dependence  of 
the  phase  boundary  movement  was  also  related  to  the  diffusion- 
controlled  growth  in  the  Johnson-Mehl-Avrami  equation  [49,69] 
and  was  also  observed  in  in  situ  measurements  of  a  single  Sn02 
nanowire  [71]. 

3.3.  Linear  sweep  voltammetry 

In  linear  sweep  voltammetry,  a  linear  scanning  voltage  Av  = 
vlsvx  was  applied  and  the  current  j(l )  was  monitored  as  a  function 
of  time.  The  boundary  condition  used  in  Eq.  (48)  was  transformed 
to  the  following: 

,Y1^  =  0)-fCl)  +  vLSVT 

JU)= - ^ - •  (o/J 

Initially,  the  concentration  X0  was  a  uniform  value  of  0.99  in  the 
solid,  and  the  simulation  was  performed  in  the  charging  process. 
The  scanning  rates  vLSV  were  0.05,  0.1,  0.2,  0.5,  and  1.  Time  was 
incremented  in  steps  of  0.02  till  X(l)  fell  below  0.01. 

Currentj(l)is  shown  as  a  function  of  voltage  ^(1,  r  =  0)  +  vLSVr 
in  Fig.  1 1  in  voltammetry  plot.  Typical  broad  peaks  were  observed 
in  the  solid  solutions,  Fig.  11(a)  forg  =  0  and  Fig.  11(b)  forg  =  -  3.  The 
peak  position  for  slow  rates  was  at  approximately  zero  voltage,  and 
it  moved  to  a  higher  voltage  at  faster  scanning  rates.  The  peak  cur¬ 
rent  also  increased  correspondingly.  The  difference  between  R  =  1 
and  R  =  0.1  was  small,  which  suggested  that  they  were  close  to  the 
reversible  (purely  diffusion-controlled)  mechanism. 
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The  voltammetry  plot  of  the  phase  transformation  material 
g  =  -  5  is  shown  in  Fig.  11(c).  For  large  R  =  10,  interface  control,  the 
plot  was  similar  to  those  of  solid  solutions.  Sudden  current  change 
was  observed  because  of  the  initiation  and  termination  of  the  phase 
transformation.  For  small  R  =  0.1,  diffusion  control,  a  sharp  current 
peak  was  observed.  Under  pure  diffusion  control,  a  small  varia¬ 
tion  of  voltage  caused  a  large  current  change  because  the  voltage 
was  flat  for  the  phase  transformation  material.  For  the  interme¬ 
diate  R  =  1  value,  a  two-step  current  decay  was  observed  at  low 
scanning  rates.  Experimental  voltammetry  plots  of  battery  mate¬ 
rials,  such  as  for  LiCo02  [28,72,73],  LiMn204  [72],  graphite  [74], 
and  LiFeP04  [75,76]  generally  have  peak  shapes  similar  to  those 
of  solid  solutions  in  Fig.  11(a)  and  (b).  However,  they  are  different 


from  those  of  phase  transformation  material  in  Fig.  1 1(c),  although 
peak  positions  in  those  experimental  plots  correspond  to  the  two- 
phase  coexistence  region  in  these  materials.  The  simulation  results 
in  Fig.  11(c)  represent  only  a  single  particle,  and  the  dispersion  of 
time  constants  in  a  composite  electrode  could  broaden  the  peaks 
in  Fig.  11(c). 

Finally,  the  scanning  rate  dependence  of  the  peak  current  data 
for  all  three  thermodynamics  is  shown  in  Fig.  11(d).  In  classical 
liquid  electrochemistry,  the  peak  current  of  a  redox  reaction  is 
proportional  to  the  square  root  of  scanning  rate  VvLSV,  known 
as  Randles-Sevcik  equation  [15,22].  This  was  derived  under  the 
condition  of  the  reversible  reaction  of  an  oxidant  and  reductant, 
semi-infinite  diffusion,  and  constant  concentration  at  infinity.  In 


R=10  R=1  R=0.1 


R=10  R=1  R=0.1 
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Fig.  11.  Linear  sweep  voltammetry  curves  under  three  different  rates  (0.05, 0.2,  and  1 )  and  interfacial  resistances  R  of  three  different  thermodynamics  (a)  g  =  0,  (b)  g  =  -  3, 
and  (c)  g=  -  5.  (d)  Peak  current  as  a  function  of  scanning  rate  for  g= 0  and  square  root  of  scanning  rate  for  g=  -  3,  -  5. 
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Fig.  11.  (  Continued ). 


the  present  work,  the  peak  current  was  approximately  propor¬ 
tional  to  vLSV  when  g  =  0  and  it  was  roughly  proportional  to  vV-w 
when  g=  -  3,  -  5.  The  clear  exception  was  the  case  of  reversible 
phase  transformation  (ft  =  0.1).  The  curves  for  ft  =  1,  10  were  sim¬ 
ilar  to  those  in  the  experimental  investigation  of  LiFeP04  [76].  In 
literatures,  both  the  linear  [75,77]  and  the  square  root  dependence 
[78,79]  on  the  scanning  rate  have  been  observed  for  LiFePOzp 

3.4.  Impedance  spectroscopy 

In  the  simulation  of  impedance  spectroscopy,  a  sinusoidal  cur¬ 
rent  j(l)  was  applied  and  voltage  response  t/t(1)  was  monitored. 
Similar  to  the  constant-current  case  in  Section  3.1,  the  additional 
interface  contribution  can  be  easily  obtained.  The  boundary  condi¬ 
tion  to  be  used  in  Eq.  (48)  became  the  following: 

K 1 )  =  J'ac  sin(<wr)  =  jAC  sin(2:r/r),  (68) 

where  the  dimensionless  frequency  was  defined  as  follows: 

/  =  /exp  to,  (69) 

where /exp  is  the  experimental  frequency  in  Hz. 

Initially,  the  concentration  X0  was  a  uniform  value  of  0.9  in 
the  solid.  The  amplitude  of  the  sine  wave  jAC  was  0.0001,  0.001, 
0.01,  and  0.05.  The  dimensionless  frequency /ranged  from  0.05  to 
1000.  For  each  frequency,  a  digital  sine  wave  with  1024  uniformly 
spaced  points  was  applied  for  one  period.  Thus,  the  time  stepping 
was  1  /(1024/)  until  it  reached  1  //.  The  translation  of  dimensionless 
frequency  /to  experimental  frequency  /exp  dependedon  t0,  as  in  Eq. 
(69).  Thus,  for  a  value  of  1 000  s  of  t0, 0.05  of/translated  to  0.05  mHz 
of  /exp,  which  was  below  the  frequency  that  is  normally  applied  in 
experiments.  A  value  of  1 000  for/translated  to  1  Hz,  which  was  low 
for  the  high  frequency  limit  in  real  experiments.  From  the  view¬ 
point  of  simulation,  this  dimensionless  high  frequency  limit  can 
be  easily  extended  to  above  1000,  but  simulation  will  require  an 
increasing  number  of  elements  to  insure  convergence.  For  example, 
5000  instead  of  500  elements  were  required  when  a  dimensionless 
frequency  of  106was  used. 

As  for  the  cases  of  constant-current  and  constant-voltage  sig¬ 
nals  in  Sections  3.1  and  3.2,  the  evolution  of  concentration  profiles 
was  presented  in  Fig.  1 2  for  a  sine  current  wave  with  dimensionless 
amplitudes  of  0.01  (left  figure)  and  0.05  (right  figure)  at  a  dimen¬ 
sionless  frequency  of  0.05,  for  the  phase  transformation  material 
with  g=-  5.  As  expected,  the  oscillation  of  current  at  the  particle 
surface  x  =  1  led  to  the  oscillation  of  surface  concentration.  For  the 
small  amplitude  of  0.01 ,  phase  separation  was  not  initiated.  For  the 
larger  amplitude  of  0.05,  the  particle  exhibited  a  two-phase  coex¬ 
istence  at  a  quarter  of  a  period  with  r  =  5,  and  the  phase  boundary 


oscillated  from  r  =  5  to  r  =  15.  The  concentration  at  the  end  of  the 
period  was  almost  the  same  as  the  initial  condition.  The  evolution 
of  concentration  profiles  in  the  solid  solutions  was  similar  to  values 
in  the  left  figure  of  Fig.  12. 

The  solution  of  Eq.  (48)  with  boundary  conditions  Eq.  (68)  pro¬ 
vided  the  time  domain  voltage  signal  A^(l)  =  r  =  0), 

which  can  be  expanded  as  a  Fourier  series  [80] 

oo 

AWl,r)=^  +  E  Ancos{ncor  +  </>n),  (70) 

n= 1 

where  ao/2  is  the  Faradaic  rectification  (DC  component),  An  and 
(pn  are  the  amplitude  and  phase  angle  of  the  nth  harmonic  term. 
The  classical  impedance  spectroscopy  only  looks  at  the  1st  or  fun¬ 
damental  harmonic  term  as  in  Eq.  (61).  These  coefficients  can  be 
obtained  from  a  Fast  Fourier  Transform  based  Fourier  (harmonic) 
analysis  detailed  in  Ref.  [80]. 

Impedance  spectra  of  the  1st  harmonic  in  the  form  of  a  Nyquist 
plot  are  shown  in  Fig.  13(a)  for  the  different  amplitudes  investi¬ 
gated,  along  with  the  analytical  results  (dashed  lines)  from  Eq.  (61 ) 
by  replacing  s  with  ico.  The  spectrum  shape  behaved  as  a  series  com¬ 
bination  of  a  resistor  and  a  capacitor  at  low  frequencies  and  has 
been  traditionally  called  the  Warburg  element  with  impermeable 
[30],  reflecting  [81],  reflective  [82],  open  [81],  and  restricted  [83] 
boundary  conditions,  as  well  as  a  finite-space  Warburg  element 
[84]  or  T  element  [85].  Numerical  results  of  the  smallest  amplitude 
of  0.0001  were  closely  matched  to  analytical  results  as  expected. 
As  the  amplitude  increased,  deviations  can  be  observed  for  both 
solid  solutions  and  phase  transformation  material.  This  nonlinear 
behavior  was  due  to  the  nonlinear  form  of  the  partial  differential 
equations  Eq.  (48),  common  to  all  three  cases.  As  g  changed  from 
0  to  -3  and  -5,  the  diffusivity  of  Eq.  (60)  decreased,  and  vertical 
tail  height  decreased  at  the  same  frequency  range.  The  nonlinear 
behavior  can  also  be  quantified  from  the  harmonic  analysis;  the 
results  are  shown  in  Fig.  13(c).  The  voltage  amplitude  ratios  of  the 
0th  and  2nd  harmonic  term  to  the  1st  harmonic  term  are  plotted. 
First,  the  contribution  of  the  0th  harmonic  term  (DC  component) 
was  always  relatively  large  but  this  information  is  rarely  studied 
theoretically  or  experimentally.  Second,  the  2nd  harmonic  term  at 
the  high  frequency  range  increased  slightly  with  increasing  ampli¬ 
tude.  In  addition,  even  at  the  smallest  amplitude  of  0.0001,  this 
contribution  was  close  to  10%,  and  a  slight  deviation  from  analyti¬ 
cal  results  (dashed  lines)  can  be  clearly  observed  in  Fig.  13(b),  with 
g=-  5  and  amplitude  0.05  as  an  example.  Third,  the  2nd  harmonic 
term  at  the  low  frequency  range  increased  greatly  with  increas¬ 
ing  amplitude  which  can  also  be  observed  by  the  clear  deviation 
in  Fig.  13(a).  Finally,  for  the  phase  transformation  material  with 
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Fig.  12.  Evolution  of  the  concentration  profile  under  a  sinusoidal  current  of  amplitude  j/\c  =  0.01  (left)  and  jAc  =  0.05  (right)  for  the  phase  transformation  material  withg=  -  5. 
Dimensionless  times  are  shown  on  each  curve.  The  simulation  was  performed  for  one  period  of  time  20.  The  curves  corresponding  to  time  20  are  almost  the  same  as  that  at 
time  0  in  both  figures.  Particle  surface  and  center  are  at  the  position  of  1  and  0,  respectively. 


g=  -  5,  a  large  sine  current  of  0.05  reached  the  unstable  region  and 
initiated  the  phase  separation  at  low  frequencies.  This  led  to  the 
largest  deviation  from  the  analytical  results,  which  can  be  seen  in 
Fig.  13(b). 

3.5.  Initial  condition  of  “X=0.5” 

In  addition  to  the  initial  concentrations  of  0.9  and  0.99  discussed 
above,  a  value  of  0.5  seemed  to  be  another  interesting  initial  con¬ 


centration.  However,  as  can  be  seen  in  Fig.  2(a),  this  concentration 
corresponds  to  a  separation  of  two  phases  for  the  phase  transforma¬ 
tion  material  withg=  -  5.  In  other  words,  a  uniform  concentration 
profile  of  0.5  is  not  a  valid  initial  condition  forg  =  -  5. 

Instead,  the  particle  was  charged  at  a  constant  current  of  0.1  for 
the  duration  of  dimensionless  time  of  4,  starting  from  a  uniform 
concentration  of  0.9.  The  delivered  capacity  will  change  the  overall 
average  concentration  to  0.5.  Then  the  current  was  turned  off  and 
the  particle  was  relaxed  for  another  time  period  of  4.  For  solid  solu- 
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Fig.  13.  (a)  Nyquist  plot  of  the  1st  order  harmonic  term  under  different  current  amplitudes  jAC  (0.0001,  0.001,  0.01,  and  0.05)  for  three  thermodynamic  parameters  g  =  0,  (b) 
g=  -  3,  and  (c)  g=  -  5.  (b)  Enlarged  Nyquist  plot  at  high  frequencies  under  an  amplitude  of  0.05  for  g=  -  5.  (c)  Voltage  amplitude  ratios  of  the  0th  and  2nd  order  harmonic 
terms  to  that  of  the  1st  order  harmonic  term  at  different  current  amplitudes  and  frequencies. 
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Fig.  14.  Evolution  of  concentration  profile  of  a  solid  solution  g  =  0  (left)  and  a  phase  transformation  material  g  =  -  5  (right),  under  a  constant  current  j(l)  =  0.1  for  a  time  of  4 
followed  by  relaxation  for  a  time  of  4.  Dimensionless  times  are  shown  on  each  curve.  Particle  surface  and  center  are  at  the  position  of  1  and  0,  respectively. 


tions,  the  concentration  profile  was  expected  to  relax  to  a  uniform 
value  of  0.5  across  the  particle.  The  case  ofg=0  is  shown  in  the  left 
figure  of  Fig.  14  with  labeled  dimensionless  times.  Similar  behavior 
was  found  for  the  case  ofg  =  -  3.  However,  forg  =  -  5  (right  figure), 
the  final  relaxed  state  of  apparent  concentration  of  0.5  exhibited  the 
coexistence  of  two  phases.  The  two  phases  had  the  same  composi¬ 
tions  of  binodal  points  B1  and  B2  in  Fig.  2(a).  These  two  different 
behaviors  can  be  understood  by  Eq.  (48).  A  zero  current  during 
relaxation  implied  a  constant  value  of  ^  across  the  particle.  For 
solid  solutions,  y  =  0  led  to  constant  values  of  concentration  every¬ 
where.  Conversely,  a  concentration  gradient  was  possible  with  a 
nonzero  value  of  y.  A  good  analog  was  the  space  charge  concept 
with  the  balance  of  concentration  and  electrical  driving  forces  [86]. 

3.5 A.  GITT 

In  the  context  of  GITT,  the  right-hand  side  of  Eq.  (63)  is  a 
constant.  This  condition  was  not  satisfied  for  an  initial  uniform  con¬ 
centration  of  0.99  shown  in  Fig.  7.  The  same  types  of  plots  are  shown 
in  Fig.  1 5  for  solid  solutions  with  a  uniform  concentration  of  0.5  and 
phase  transformation  material  with  initial  two-phase  coexistence. 
The  dashed  lines  are  the  calculated  results  from  Eq.  (63).  A  plateau 
was  obtained  up  to  an  approximate  time  of  0.5  forg  =  0  and  1  for 
g=  -  3.  The  data  with  g=  -  5  were  close  to  constant  and  oscillated 
around  the  dashed  line. 

The  difference  between  Figs.  7  and  15  can  be  understood  by  the 
assumption  used  to  obtain  Eq.  (63),  i.e.  the  application  of  current 


step  leading  to  a  small  value  of  A\j/.  For  the  initial  concentration  of 
0.99,  the  equilibrium  curve  was  steeper  compared  to  those  of  the 
initial  concentration  of  0.5,  Fig.  2(a).  Thus,  it  was  more  difficult  to 
meet  the  small-signal  approximation  for  concentration  of  0.99. 

3.5.2.  PITT 

For  the  initial  concentration  profile  in  Fig.  14  of  the  phase  trans¬ 
formation  material,  current  decay  following  the  voltage  step  is 
shown  in  Fig.  16.  Compared  to  the  initial  condition  with  a  single¬ 
phase,  Fig.  10(e),  the  first  stage  single-phase  decay  and  current 
spike  disappeared  because  the  particle  was  initially  at  the  two- 
phase  coexistence. 

3.5.3.  EIS 

For  solid  solutions  with  g  =  0,  -3,  simulation  results  for  the 
initial  uniform  concentration  of  0.5  (Fig.  17(a)),  exhibited  similar 
shapes  compared  to  those  for  the  concentration  of  0.9,  Fig.  13(a). 
Again,  dashed  lines  represented  results  predicted  from  the  analyt¬ 
ical  expression  of  Eq.  (61).  Impedance  spectra  appeared  linear  up 
to  an  amplitude  of  0.01.  At  higher  amplitude,  e.g.  0.05,  impedance 
values  became  larger  as  opposed  to  smaller  in  Fig.  13(a).  Harmonic 
analysis,  Fig.  17(b),  revealed  similar  DC  and  2nd  harmonic  terms 
compared  to  those  in  Fig.  13(c). 

However,  in  Fig.  17(a),  the  impedance  spectra  of  the  phase 
transformation  material  g=-  5  have  different  shapes  from  those 
for  the  concentration  of  0.9,  Fig.  13(a).  This  shape  behaved  as  a 


Fig.  15.  Short-time  characteristic  curves  according  to  Eq.  (63)  under  three  different  currents  for  three  different  thermodynamic  parameters  g= 0,  g=  -  3,  and  g=  -  5.  Dashed 
lines  represent  the  plateaus  expected  from  the  small-signal  assumption.  The  initial  condition  is  the  relaxed  state  in  Fig.  14,  i.e.  “X=0.5”. 
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Fig.  16.  Current  response  for  phase  transformation  material  withg  =  -  5  under  three  different  voltage  step  and  interfacial  resistance  values.  The  initial  condition  is  the  relaxed 
state  in  Fig.  14,  i.e.  “X=0.5”. 


resistor  instead  of  a  capacitor  at  low  frequencies  and  has  been 
traditionally  called  the  Warburg  element  with  Nernstian  [30], 
absorbing  [81],  transmissive  [82],  short  [81],  and  bounded  [83] 
boundary  conditions,  and  finite-length  Warburg  element  [84]  or 
O  element  [85].  The  dashed  line  corresponds  to  a  reflective  War¬ 
burg  element  predicted  according  to  Eq.  (61)  from  a  uniform 
initial  concentration  of  the  binodal  point  Bl.  The  outer  half  of  the 
particle  had  this  concentration  of  Bl  as  can  be  seen  in  Fig.  14. 
Similar  experimental  impedance  spectra  compared  to  those  in 
Fig.  17(a)  were  observed  for  the  50%  charged  LiFeP04  electrode 
[87]. 


The  analytical  result  that  led  to  the  reflective  Warburg  element 
was  obtained  under  the  assumption  of  a  constant  initial  concen¬ 
tration  and  a  small-signal  perturbation.  Thus,  it  was  natural  to 
compare  concentration  profiles  obtained  from  the  current  initial 
condition  and  those  from  the  uniform  initial  concentration  of  0.9. 
Fig.  18  shows  the  concentration  evolution  under  a  sine  wave  of 
amplitudes  0.01  (left)  and  0.05  (right)  at  a  frequency  of  0.05.  The 
initial  condition  r  =  0  was  the  relaxed  state  in  Fig.  14.  It  was  appar¬ 
ent  that  neither  of  the  figures  in  Fig.  18  satisfied  the  constant  initial 
concentration  condition,  although  the  signal  on  the  left  figure  can 
be  considered  small  amplitude.  Oscillation  of  surface  concentration 
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Fig.  17.  (a)  Nyquist  plot  of  the  1st  order  harmonic  term  under  different  current  amplitudes  for  three  thermodynamic  parameters  g=0,  (b)  g=  -  3,  and  (c)  g=  -  5.  The  initial 
condition  is  the  relaxed  state  in  Fig.  14,  i.e.  “X=0.5”.  (b)  Voltage  amplitude  ratios  of  the  0th  and  2nd  order  harmonic  terms  to  that  of  the  1st  order  harmonic  term  at  different 
current  amplitudes  and  different  frequencies. 
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Fig.  18.  Evolution  of  the  concentration  profile  of  a  solid  solution  g  =  0  and  a  phase  transformation  material  g=  -  5,  under  sine  waves  of  amplitudes  0.01  (left)  and  0.05  (right) 
at  a  frequency  of  0.05.  The  initial  condition  is  the  relaxed  state  in  Fig.  14,  i.e.  “X=0.5”.  Dimensionless  times  are  shown  on  each  curve.  The  simulation  was  performed  for 
one  period  of  time  20.  The  curve  corresponding  to  time  20  is  almost  the  same  as  that  at  time  0  for  the  left  figure.  Particle  surface  and  center  are  at  the  position  of  1  and  0, 
respectively. 


was  accompanied  by  the  “in-phase”  oscillation  of  phase  bound¬ 
ary,  and  this  was  the  reason  that  a  transmissive  Warburg  element 
was  observed  in  Fig.  17(a)  (g=  -  5).  Harmonic  analysis  of  Fig.  17(b) 
reveals  larger  harmonic  terms  with  higher  input  amplitudes  at  low 
frequencies.  Finally,  it  was  interesting  to  note  that  a  large  sine  cur¬ 
rent  could  induce  a  peculiar  phase  I|II|I  sandwich  structure,  e.g. 
between  time  15  and  20  on  the  right  figure  of  Fig.  18.  It  is  possible 
that  this  structure  can  also  be  obtained  by  applying  any  large  ampli¬ 
tude  discharging  signal  from  the  initial  condition  in  Fig.  14,  although 
this  has  been  rarely  discussed  theoretically  or  experimentally. 


forms  for  the  voltage  step  with  R  =  10  and  the  linear  sweep  voltage 
with  R=  1  are  shown  in  Fig.  19.  In  both  figures,  the  distinction 
of  current  decay  due  to  phase  boundary  movement  and  the  fol¬ 
lowing  single-phase  diffusion  was  subtler  for  spherical  symmetry 
compared  to  planar  symmetry.  Another  difference  was  the  time 
dependence  of  the  phase  boundary  position  in  different  symme¬ 
tries.  The  phase  boundary  position  to  the  second  (spherical),  one 
and  half  (cylindrical),  and  first  (planar)  power  had  a  square  root 
and  linear  time  dependence  under  the  voltage  and  current  step, 
respectively. 


3.6.  Symmetry  effect 

All  the  above  numerical  results  were  obtained  for  the  planar 
symmetry  m  =  0  to  facilitate  a  comparison  with  the  well-known 
relations  such  as  the  Cottrell,  Randles-Sevcik  equation,  etc.  Small- 
signal  solutions  with  short-time  and  long-time  approximations 
have  been  determined  for  spherical,  cylindrical,  and  planar  symme¬ 
try,  under  the  application  of  a  constant-current,  constant-voltage, 
and  sinusoidal  current/voltage  in  Ref.  [12].  The  numerical  simula¬ 
tion  with  arbitrary  magnitude  signal  suggested  that  many  features 
shown  in  previous  sections  were  common  to  all  three  symmetries 
(not  shown).  One  difference  worth  pointing  out  was  that  the  two- 
stage  current-decaying  behavior  found  for  the  constant-voltage 
and  linear  sweep  voltage  signals  (Figs.  10(e)  and  11(c)),  was  less 
distinct  for  the  cylindrical  and  spherical  symmetry.  Current  wave- 


3.7.  Applicability  of  present  work  to  real  battery  materials 

This  simulation  study  provided  herein  presents  a  framework 
that  employs  a  regular  solution  model  (thermodynamic)  and  a 
generalized  Poisson-Nernst-Planck  equation  set  with  a  resistive 
boundary  condition  (kinetic)  on  a  single  particle  (microstructure). 
The  parameter  space  was  [c° ,  r0,  /x+ ,  g,  y,  ,  R].  Variables  c+  and 
r0  were  the  maximum  concentration  and  particle  size  along  the 
symmetry  direction,  respectively.  Variables  /x°  and  g  were  related 
to  the  voltage  plateau  and  miscibility  gap  [47]  of  the  voltage  curve. 
Variable  y  represented  the  interfacial  energy  between  coexisting 
solid  phases  during  the  phase  transformation.  Diffusivity  was 
the  primary  kinetic  parameter  which  was  also  related  to  the  ionic 
conductivity  or  mobility  through  Eqs.  (27)  and  (28).  Finally,  R  rep¬ 
resented  the  resistance  across  the  solid  [liquid  interface. 


Fig.  19.  (a)  Current  response  for  the  spherical  phase  transformation  material  with  g=- 5  under  three  different  voltage  steps  Av  (0.2,  0.5,  and  1)  and  R=  10.  The  initial 
condition  is  a  uniform  concentration  of  0.9.  (b)  Linear  sweep  voltammetry  curves  for  spherical  phase  transformation  material  g=  -  5  under  three  different  rates  (0.05, 
0.2,  and  1)  and  R=  1.  The  initial  condition  is  a  uniform  concentration  of  0.99. 
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One  or  more  components  of  the  thermodynamics,  kinetics, 
and  microstructure  may  have  to  be  improved  to  allow  more 
quantitative  comparison  with  experimental  data,  usually  from 
composite  electrodes.  For  example,  the  regular  solution  model 
can  only  roughly  approximate  the  phase  transformation  materials 
such  as  LixFeP04  [63,88]  and  Li4+3XTi50i2  [89],  with  pseudo  anti- 
symmetrical  voltage  curves.  In  the  previous  phase  field  study  of 
LixFeP04,  g  values  of  -4.7  [39]  and  -10  [38,48]  were  used.  The 
interfacial  energy  between  coexisting  phase  is  normally  difficult  to 
measure  experimentally,  and  estimates  of  values  related  to  inter¬ 
facial  energy  parameter  y  were  also  provided  in  references  for 
LixFeP04  [38,39,48].  On  the  other  hand,  the  regular  solution  model 
generally  cannot  be  applied  to  other  common  materials,  such  as 
LixCo02  [22],  LixC6  [90,91],  or  LixMn204  [92,93],  which  have  more 
complex  voltage  dependences.  These  materials  have  more  than  one 
two-phase  coexistence  regions.  A  discussion  on  the  correlation  of 
other  thermodynamic  models  and  properties  with  practical  battery 
materials  will  appear  in  another  manuscript. 

The  present  simulation  can  easily  be  extended  to  more  than 
one  dimension  by  using  the  original  equations  with  gradients  and 
Laplacians,  Eqs.  (5)-(8).  In  addition,  the  single  particle  considera¬ 
tion  can  easily  be  extended  to  a  collection  of  particles,  where  either 
pseudo  2D  symmetry  from  volume  averaging  [13,51,52],  or  the 
brute-force  2D  [94]  or  3D  [95]  calculations  can  be  performed. 

4.  Conclusions 

Numerical  modeling  of  a  single  intercalation  battery  particle 
was  performed  in  this  study  for  four  different  electrical  signals, 
such  as  the  current  step,  voltage  step,  linear  sweep  voltage,  and  the 
sinusoidal  voltage.  Techniques  that  utilize  these  signals  are  gen¬ 
erally  known  as  constant-current  charging/GITT,  constant-voltage 
charging/PITT,  CV,  and  EIS  in  the  battery  field.  Thermodynamic  and 
kinetic  properties  of  interest  included  a  solid  solution  with  con¬ 
stant  diffusivity,  a  solid  solution  with  variable  diffusivity,  and  a 
phase  transformation  material.  In  addition,  the  interfacial  behav¬ 
ior  of  particle  and  liquid  electrolyte  was  incorporated  as  a  resistor. 
It  was  found  that  the  known  small-signal  solutions  of  GITT,  PITT, 
and  EIS  were  valid  only  under  the  conditions  of  a  uniform  initial 
concentration  profile,  a  small-signal  perturbation,  and  no  interfa¬ 
cial  effects.  Some  interesting  phenomena  were  observed  from  the 
simulation  of  the  phase  transformation  material.  Typical  progres¬ 
sion  for  the  current/voltage  response  in  the  phase  transformation 
exhibited  a  single-phase  diffusion,  the  initiation  of  phase  separa¬ 
tion,  movement  of  the  phase  boundary,  and  another  single-phase 
diffusion.  This  has  important  implications  in  the  interpretation  of 
experimental  data  involving  phase  transformation  battery  materi¬ 
als. 

Acknowledgement 

W.L.  would  like  to  thank  Dr.  Francesco  Ciucci  for  his  helpful 
comments  on  the  manuscript  and  Michigan  State  University  for 
providing  the  start-up  package. 

References 

[  1  ]  R.P.  Buck,  J.  Membr.  Sci.  1 7  ( 1 984)  1  -62. 

[2]  M.Z.  Bazant,  K.  Thornton,  A.  Ajdari,  Phys.  Rev.  E  70  (2004)  021506. 

[3]  T.  Grasser,  T.W.  Tang,  H.  Kosina,  S.  Selberherr,  Proc.  IEEE  91  (2003)  251-274. 

[4]  R.S.  Eisenberg,  J.  Membr.  Biol.  150  (1996)  1-25. 

[5]  J.  Crank,  The  Mathematics  of  Diffusion,  2nd  ed.,  Oxford  University  Press,  New 
York,  NY,  1980. 

[6]  T.R.  Brumleve,  R.P.  Buck,  J.  Electroanal.  Chem.  126  (1981)  73-104. 

[7]  R.P.  Buck,  C.  Mundt,  Electrochim.  Acta  44  (1999)  1999-2018. 

[8]  J.  Jamnik,  J.  Maier,  Phys.  Chem.  Chem.  Phys.  3  (2001)  1668-1678. 

[9]  W.  Lai,  S.M.  Haile,  J.  Am.  Ceram.  Soc.  88  (2005)  2979-2997. 

[10]  W.  Lai,  S.M.  Haile,  Phys.  Chem.  Chem.  Phys.  10  (2008)  865-883. 


[11]  W.  Lai,  F.  Ciucci,  Electrochim.  Acta  56  (2010)  531-542. 

[12]  W.  Lai,  F.  Ciucci,  J.  Electrochem.  Soc.  158  (2011)  A115-A121. 

[13]  W.  Lai,  F.  Ciucci,  Electrochim.  Acta.  doi:10.1016/j.electacta.2011.01.012. 

[14]  A.J.  Bard,  L.R.  Faulkner,  Electrochemical  Methods:  Fundamentals  and  Applica¬ 
tions,  2nd  ed.,  John  Wiley  &  Sons  Inc.,  New  York,  NY,  2000. 

[15]  A.J.  Bard,  G.  Inzelt,  F.  Scholz  (Eds.),  Electrochemical  Dictionary,  Springer,  Berlin, 
Germany,  2008. 

[16]  D.  Linden,  T.B.  Reddy  (Eds.),  Handbook  of  Batteries,  McGraw-Hill  Companies 
Inc.,  New  York,  NY,  2001. 

[17]  W.  Weppner,  R.A.  Huggins,  J.  Electrochem.  Soc.  124  (1977)  1569-1578. 

[18]  C.J.  Wen,  B.A.  Boukamp,  R.A.  Huggins,  W.  Weppner,  J.  Electrochem.  Soc.  126 
(1979)2258-2266. 

[19]  E.  Barsoukov,  J.R.  Macdonald  (Eds.),  Impedance  Spectroscopy:  Theory,  Experi¬ 
ment,  and  Applications,  John  Wiley  &  Sons  Inc.,  Hoboken,  NJ,  2005. 

[20]  F.  Kremer,  A.  Schonhals  (Eds.),  Broadband  Dielectric  Spectroscopy,  Springer, 
Berlin,  Germany,  2002. 

[21  ]  M.E.  Orazem,  B.  Tribollet,  Electrochemical  Impedance  Spectroscopy,  John  Wiley 
&  Sons  Inc.,  Hoboken,  NJ,  2008. 

[22]  M.D.  Levi,  G.  Salitra,  B.  Markovsky,  H.  Teller,  D.  Aurbach,  U.  Heider,  L.  Heider,  J. 
Electrochem.  Soc.  146  (1999)  1279-1289. 

[23]  S.  Atlung,  K.  West,  T.  Jacobsen,  J.  Electrochem.  Soc.  126  (1979)  1311-1321. 

[24]  E.  Deiss,  Electrochim.  Acta  50  (2005)  2927-2932. 

[25]  V.R.  Subramanian,  J.A.  Ritter,  R.E.  White,  J.  Electrochem.  Soc.  148  (2001) 
E444-E449. 

[26]  C.  Montella,  J.  Electroanal.  Chem.  518  (2002)  61-83. 

[27]  E.  Deiss,  Electrochim.  Acta  47  (2002)  4027-4034. 

[28]  H.C.  Shin,  S.I.  Pyun,  Electrochim.  Acta  46  (2001)  2477-2485. 

[29]  C.  Montella,  J.  Electroanal.  Chem.  614  (2008)  121-130. 

[30]  T.  Jacobsen,  K.  West,  Electrochim.  Acta  40  (1995)  255-262. 

[31  j  S.R.  De  Groot,  P.  Mazur,  Non-Equilibrium  Thermodynamics,  Dover  Publications, 
Mineola,  NY,  1984. 

[32]  F.  Ciucci,  W.  Lai,  Transport  Porous  Media.  doi:10.1007/sl  1242-01  lr-r9738-5. 

[33]  P.  Atkins,  J.  De  Paula,  Physical  Chemistry,  Oxford  University  Press,  Oxford,  UK, 
2006. 

[34]  S.  Stolen,  T.  Grande,  N.L.  Allan,  Chemical  Thermodynamics  of  Materials,  John 
Wiley  &  Sons,  West  Sussex,  UK,  2004. 

[35]  W.R.  McKinnon,  R.R.  Haering,  in:  R.E.  White,  J.O.M.  Bockris,  B.E.  Conway  (Eds.), 
Modern  Aspects  of  Electrochemistry,  vol.  15,  Plenum  Press,  New  York,  NY,  1983, 
pp.  235-304. 

[36]  S.T.  Coleman,  W.R.  McKinnon,  J.R.  Dahn,  Phys.  Rev.  B  29  (1984)  4147-4149. 

[37]  M.D.  Levi,  D.  Aurbach,  Electrochim.  Acta  45  (1999)  167-185. 

[38]  G.K.  Singh,  G.  Ceder,  M.Z.  Bazant,  Electrochim.  Acta  53  (2008)  7599-7613. 

[39]  B.C.  Han,  A.  Van  der  Ven,  D.  Morgan,  G.  Ceder,  Electrochim.  Acta  49  (2004) 
4691-4699. 

[40]  H.B.  Callen,  Thermodynamics  and  an  Introduction  to  Thermostatistics,  John 
Wiley  &  Sons,  1985. 

[41  ]  K.  Thornton,  J.  Agren,  P.W.  Voorhees,  Acta  Mater.  51  (2003)  5675-5710. 

[42]  L.Q,  Chen,  Ann.  Rev.  Mater.  Res.  32  (2002)  113-140. 

[43]  M.  Tang,  W.C.  Carter,  R.M.  Cannon,  Phys.  Rev.  B  73  (2006)  024102. 

[44]  W.J.  Boettinger,  J.A.  Warren,  C.  Beckermann,  A.  Karma,  Ann.  Rev.  Mater.  Res.  32 
(2002)163-194. 

[45]  J.E.  Guyer,  W.J.  Boettinger,  J.A.  Warren,  G.B.  McFadden,  Phys.  Rev.  E  69  (2004) 
021603. 

[46]  J.E.  Guyer,  W.J.  Boettinger,  J.A.  Warren,  G.B.  McFadden,  Phys.  Rev.  E  69  (2004) 
021604. 

[47]  D.  Burch,  M.Z.  Bazant,  Nano  Lett.  9  (2009)  3795-3800. 

[48]  M.  Tang,  H.Y.  Huang,  N.  Meethong,  Y.H.  Kao,  W.C.  Carter,  Y.M.  Chiang,  Chem. 
Mater.  21  (2009)  1557-1571. 

[49]  M.  Tang,  W.C.  Carter,  Y.M.  Chiang,  Ann.  Rev.  Mater.  Res.  40  (2010)  501-529. 

[50]  J.W.  Cahn,  J.E.  Hilliard,  J.  Chem.  Phys.  28  (1958)  258-267. 

[51]  M.  Doyle,  T.F.  Fuller,  J.  Newman,  J.  Electrochem.  Soc.  140  (1993)  1526-1533. 

[52]  T.F.  Fuller,  M.  Doyle,  J.  Newman,  J.  Electrochem.  Soc.  141  (1994)  1-10. 

[53]  J.  Newman,  K.E.  Thomas-Alyea,  Electrochemical  Systems,  John  Wiley  &  Sons 
Inc.,  Hoboken,  NJ,  2004. 

[54]  I.  Riess,  J.  Maier,  J.  Electrochem.  Soc.  156  (2009)  P7-P20. 

[55]  H.C.  Chang,  G.  Jaffe,  J.  Chem.  Phys.  20  (1952)  1071-1077. 

[56]  A.V.  Virkar,  J.  Power  Sources  147  (2005)  8-31. 

[57]  A.V.  Virkar,  J.  Power  Sources  194  (2009)  753-762. 

[58]  H.C.  Shin,  S.I.  Pyun,  Electrochim.  Acta  45  (1999)  489-501. 

[59]  COMSOL  Multiphysics,  1998-2011  COMSOLAB. 

[60]  J.  Molenda,  A.  Stoklosa,  T.  Bak,  Solid  State  Ionics  36  (1989)  53-58. 

[61  ]  W.  Lai,  C.K.  Erdonmez,  T.F.  Marinis,  C.K.  Bjune,  N.J.  Dudney,  F.  Xu,  R.  Wartena, 
Y.M.  Chiang,  Adv.  Mater.  22  (2010)  E139-E144. 

[62]  Z.H.  Chen,  J.R.  Dahn,  Electrochim.  Acta  49  (2004)  1079-1090. 

[63]  N.  Meethong,  H.Y.S.  Huang,  W.C.  Carter,  Y.M.  Chiang,  Electrochem.  Solid-State 
Lett.  10  (2007)  A134-A138. 

[64]  R.W.  Balluffi,  S.M.  Allen,  W.C.  Carter,  Kinetics  of  Materials,  John  Wiley  and  Sons, 
Hoboken,  NJ,  2005,  pp.  501-531. 

[65]  D.  Aurbach,  M.D.  Levi,  E.  Levi,  Solid  State  Ionics  179  (2008)  742-751. 

[66]  N.  Meethong,  Y.H.  Kao,  W.C.  Carter,  Y.M.  Chiang,  Chem.  Mater.  22  (2010) 
1088-1097. 

[67]  A.  Funabiki,  M.  Inaba,  T.  Abe,  Z.  Ogumi,  J.  Electrochem.  Soc.  146  (1999) 
2443-2448. 

[68]  M.D.  Levi,  E.  Markevich,  D.  Aurbach,  Electrochim.  Acta  51  (2005)  98-110. 

[69]  J.L.  Allen,  T.R.  Jow,  J.  Wolfenstine,  Chem.  Mater.  19  (2007)  2108-2111. 

[70]  H.C.  Shin,  S.I.  Pyun,  S.W.  Kim,  M.H.  Lee,  Electrochim.  Acta  46  (2001)  897-906. 


W.  Lai  /  Journal  of  Power  Sources  196(2011)  6534-6553 


6553 


[71]  J.Y.  Huang,  L.  Zhong,  C.M.  Wang,  J.P.  Sullivan,  W.  Xu,  L.Q,  Zhang,  S.X.  Mao,  N.S. 
Hudak,  X.H.  Liu,  A.  Subramanian,  H.Y.  Fan,  L.A.  Qi,  A.  Kushima,  J.  Li,  Science  330 
(2010)1515-1520. 

[72]  D.A.  Totir,  B.D.  Cahan,  D.A.  Scherson,  Electrochim.  Acta  45  (1999)  161-166. 

[73]  J.  Barker,  R.  Pynenburg,  R.  Koksbang,  M.Y.  Saidi,  Electrochim.  Acta  41  (1996) 
2481-2488. 

[74]  M.D.  Levi,  D.  Aurbach,  J.  Electroanal.  Chem.  421  (1997)  79-88. 

[75]  F.  Sauvage,  E.  Baudrin,  L.  Gengembre.J.M.Tarascon,  Solid  State  Ionics  176  (2005) 
1869-1876. 

[76]  D.Y.W.  Yu,  C.  Fietzek,  W.  Weydanz,  K.  Donoue,  T.  Inoue,  H.  Kurokawa,  S.  Fujitani, 
J.  Electrochem.  Soc.  1 54  (2007)  A253-A257. 

[77]  S.  Franger,  C.  Bourbon,  F.  Le  Cras, J.  Electrochem.  Soc.  151  (2004)  A1024-A1 027. 

[78]  S.W.  Song,  R.P.  Reade,  R.  Kostecki,  K.A.  Striebel,  J.  Electrochem.  Soc.  153  (2006) 
A12-A19. 

[79]  M.  Takahashi,  S.  Tobishima,  K.  Takei,  Y.  Sakurai,  Solid  State  Ionics  148  (2002) 
283-289. 

[80]  W.  Lai,  Electrochim.  Acta  55  (2010)  5511-5518. 

[81  ]  J.  Bisquert,  J.  Phys.  Chem.  B  106  (2002)  325-333. 


[82]  A.  Lasia,  in:  B.E.  Conway,  J.  Bockris,  R.E.  White  (Eds.),  Modern  Aspects  of  Electro¬ 
chemistry,  vol.  32,  Kluwer  Academic/Plenum  Publishers,  New  York,  NY,  1999, 
pp.  143-248. 

[83]  J.P.  Diard,  B.  Le  Gorrec,  C.  Montella,  J.  Electroanal.  Chem.  471  (1999)  126-131. 

[84]  M.D.  Levi,  C.  Wang,  D.  Aurbach,  J.  Electroanal.  Chem.  561  (2004)  1-11. 

[85]  V.  Freger,  S.  Bason,  J.  Membr.  Sci.  302  (2007)  1-9. 

[86]  S.  Kim,  J.  Maier,  J.  Electrochem.  Soc.  149  (2002)  J73-J83. 

[87]  M.  Gaberscek,  R.  Dominko,  J.  Jamnik,  J.  Power  Sources  174  (2007)  944-948. 

[88]  A.K.  Padhi,  K.S.  Nanjundaswamy,  J.B.  Goodenough,  J.  Electrochem.  Soc.  144 
(1997)1188-1194. 

[89]  T.  Ohzuku,  A.  Ueda,  N.  Yamamoto,  J.  Electrochem.  Soc.  142  (1995)  1431-1435. 

[90]  J.R.  Dahn,  Phys.  Rev.  B  44  (1991)  9170-9177. 

[91]  T.  Ohzuku,  Y.  Iwakoshi,  K.  Sawai,  J.  Electrochem.  Soc.  140  (1993)  2490-2498. 

[92]  T.  Ohzuku,  M.  Kitagawa,  T.  Hirai,  J.  Electrochem.  Soc.  137  (1990)  769-775. 

[93]  M.M.  Thackeray,  Prog.  Solid  State  Chem.  25  (1997)  1-71. 

[94]  M.  Smith,  R.E.  Garcia,  Q.C.  Horn,  J.  Electrochem.  Soc.  156  (2009)  A896-A904. 

[95]  C.W.  Wang,  A.M.  Sastry,  J.  Electrochem.  Soc.  154  (2007)  A1035-A1047. 


